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Preface 


In nature, it is possible to observe a cooperative behaviour in all animals, since, according 
to Charles Darwin's theory, every being, from ants to human beings, form groups in which 
most individuals work for the common good. However, although study of dozens of social 
species has been done for a century, details of how and why cooperation evolved remain 
to be worked out. Actually, cooperative behaviour has been studied from different points 
of view. For instance evolutionary biologists and animal behaviour researchers look for the 
genetic basis and molecular drivers of this kind of behaviours, as well as the physiological, 
environmental, and behavioural impetus for sociality; while neuroscientists discover key 
correlations between brain chemicals and social strategies. From a more mathematical point 
of view, economics have developed a modelling approach, based on game theory, to quantify 
cooperation and predict behavioural outcomes under different circumstances. Although 
game theory has helped to reveal an apparently innate desire for fairness, developed models 
are still imperfect. Furthermore, social insect behaviour, from a biological point of view, might 
be emulated by a micro-robot colony and, in that way, analysis of a tremendous amount of 
insect trajectories and manual event counting is replaced by tracking several miniature robots 
on a desktop table. 

Swarm robotics is a new approach that emerged on the field of artificial swarm intelligence, as 
well as the biological studies of insects (i.e. ants and other fields in nature) which coordinate 
their actions to accomplish tasks that are beyond the capabilities of a single individual. In 
particular, swarm robotics is focused on the coordination of decentralised, self-organised 
multi-robot systems in order to describe such a collective behaviour as a consequence of local 
interactions with one another and with their environment. 

Research in swarm robotics involves from robot design to their controlling behaviours, by 
including tracking techniques for systematically studying swarm-behaviour. Moreover, 
swarm robotic-based techniques can be used in a number of applications. This is, for instance, 
the case of the Particle Swarm Optimization (PSO) which is a direct search method, based 
on swarm concepts, that models and predicts social behaviour in the presence of objectives. 
In this case, the swarm under study is typically modelled by particles in multidimensional 
space that have two essential reasoning capabilities: their memory of their own best position 
and the knowledge of the global or their neighbourhood's best, such that swarm members 
communicate good positions to each other and adjust their own position and velocity based 
on those good positions in order to obtain the best problem solution. 

Different challenges have to be solved in the field of swarm robotics. This book is focused 
on real practical applications by analyzing how individual robotic agents should behave in a 
robotic swarm in order to achieve a specific goal such as target localization or path planning. 
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In this context, the first paper, by Hereford and Siebold, concentrates on looking for a target 
in a room. They describe, on the one hand, the way a PSO algorithm, based on bird flocking, 
may be embedded into a robot swarm; and, on the other, the implementation of a four-step 
trophallactic behaviour of social insects in a robotic platform by making sensor measurements 
instead of exchanging information when two or more particles are in contact. Different 
software and hardware tests were developed to evaluate both search strategy performances. 

Another issue which may be solved by PSO methods is the robotic cell problem, where each 
integrating machine could be identified as a member of a swarm. In this context, Kamalabadi 
et al. present a hybrid PSO algorithm to find a schedule robot movement to minimize cycle 
time when multiple-type parts three-machine robotic cells are considered. Its performance 
has been compared with three well-known metaheuristic algorithms: Genetic Algorithm 
(GA), Basic Algorithm (PSO- 

I) and Constriction Algorithm (PSO-II), by succeeding in the most problems, especially for 
large-sized ones. 

The next two papers have focused on the problem of path planning for mobile robots. Firstly, 
Curkovic et al. introduce a way to solve the navigation problem for a robot in a workspace 
containing differently shaped and distributed by means of a simplification of Honey Bees 
Mating Algorithm. Moreover, a plan optimization technique that results in a minimization 
of the required time or the travelled distance is proposed. Again, method performance is 
successfully evaluated with respect to the Genetic Algorithm. Secondly, Xue et al. apply 
PSO-type control for real-time path planning on a typical swarm of wheeled mobile robots 
in an unstructured environment. Furthermore, an overview of a system modelling at both 
individual and swarm levels as well as a fusion-framework is presented. Their study was 
tested through virtual signal generators and simulations about swarm-component measuring 
and fusing. 

Another application of the PSO techniques is the design of an infinite impulse response (HR) 
digital filter of robot force/ position sensors. Zhang proposes an HR filter design from the 
knowledge of the structure of a filter and master of an intelligent optimization algorithm. 
The PSO algorithm is then used to optimize parameter values. Newly, simulation is used to 
validate the developed technique. 

Finally, it is essential to systematically study and test swarm-behaviour by analyzing what 
each swarm member is doing as well as where and when it acts. For that purpose, Martinez 
and del Pobil has developed a visual application that robustly identifies and tracks all robotic 
swarm members. Different situations and visual systems were studied to achieve that goal. 
Experimental results on a real system are also presented. 

This book has only provided a partial picture of the field of swarm robotics by focusing on 
practical applications. The global assessment of the contributions contained in this book is 
reasonably positive since they highlighted that it is necessary to adapt and remodel biological 
strategies to cope with the added complexity and problems that arise when robot individuals 
are considered. 


Ester Martinez Martin 
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Bio-inspired search strategies for robot swarms 

James M. Hereford and Michael A. Siebold 

Murray State University , Murray , KY 

USA 


1. Introduction 

Our goal is as follows: build a suite/ swarm of (very) small robots that can search a room for 
a "target". We envision that the robots will be about the size of a quarter dollar, or smaller, 
and have a sensor or sensors that "sniff" out the desired target. For example, the target 
could be a bomb and the robot sensors would be a chemical detectors that can distinguish 
the bomb from its surroundings. Or the target could be a radiation leak and the sensors 
would be radiation detectors. In each search scenario, we assume that the target gives off a 
diffuse residue that can be detected with a sensor. 

It is not very efficient to have the suite of robots looking randomly around the room hoping 
to "get lucky" and find the target. There needs to be some way to coordinate the 
movements of the many robots. There needs to be an algorithm that can guide the robots 
toward promising regions to search while not getting distracted by local variations. The 
search algorithm must have the following constraints: 

• The search algorithm should be distributed among the many robots. If the algorithm is 
located in one robot, then the system will fail if that robot fails. 

• The search algorithm should be computationally simple. The processor on each bot is 
small, has limited memory, and there is a limited power source (a battery) so the 
processor needs to be power efficient. Therefore, the processor will be a simple 
processor. 

• The algorithm needs to be scalable from one robot up to 10's, 100's, even 1000's of 
robots. The upper limit on the number of robots will be set by the communication links 
among the robots; there needs to be a way to share information among the robots 
without requiring lots of communication traffic. 

• The search algorithm must allow for contiguous movement of the robots. 

This chapter will describe two search strategies for robot swarms that are based on 
biological systems. The first search strategy is based on the flocking behavior of birds and 
fishes. This flocking behavior is the inspiration behind the Particle Swarm Optimization 
(PSO) algorithm that has been used in software for many types of optimization problems. In 
the PSO algorithm the potential solutions, called particles, "fly" through the problem space 
by following some simple rules. All of the particles have a fitness value based on the value 
or measurement at the particle's position and have velocities which direct the flight of the 
particles. The velocity of each particle is updated based on the particle's current velocity as 
well as the best fitness of any particle in the group. 
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We describe how the PSO algorithm can be embedded into a robot swarm by letting each 
bot's behavior be like a particle in the PSO. We call the algorithm the physically embedded 
PSO (pePSO). The bots swarm throughout the search space and take measurements. Over 
time, they cluster near the peak(s) or targets. We show through both 2D simulation results 
and robot hardware results that the pePSO effectively finds the targets with a minimum 
number of bot-bot communications. 

The second search strategy is based on the trophallactic behavior of social insects. 
Trophallaxis is the exchange of fluid by direct mouth-to-mouth contact. This phenomenon is 
observed in ants, bees, wasps and even dogs and birds. In our trophallaxis-based algorithm, 
the bots do not actually exchange information but instead make sensor measurements when 
two or more bots/ particles are "in contact". The bots remain stationary for a certain time 
that is proportional to the measurement value. Bots thus cluster in areas of the search space 
that have high fitness/ measurement values. 

This new trophallaxis-based search algorithm has several advantages over other swarm- 
based search techniques. First, no bot-bot communication is required. Thus, there is no 
concern with communication radius, protocol, or bandwidth. Second, the bots do not have 
to know their position. During the search, the hot/ particle moves randomly except when it 
collides and stops, takes a measurement, and waits. At the end of the search, the cluster 
locations can be determined from a remote camera, special-purpose robot, or human 
canvassing. 

This paper is organized as follows: section 2 gives background on the Particle Swarm 
Optimization algorithm and its use in robot swarms and section 3 gives results from 
simulations and hardware results of embedding the PSO into a robot swarm. Section 4 
discusses the trophallaxis- based search algorithm and section 5 gives simulation results 
using the trophallaxis algorithm. In section 6 we give our conclusions. 

2. Particle Swarm Optimization and robots 

In Particle Swarm Optimization (PSO) (Eberhart & Kennedy, 1995; Clerc &Kennedy, 2002; 
Eberhart & Shi, 2004), the potential solutions, called particles, "fly" through the problem 
space by following some simple rules. All of the particles have fitness values based on their 
position and have velocities which direct the flight of the particles. PSO is initialized with a 
group of random particles (solutions), and then searches for optima by updating 
generations. In every iteration, each particle is updated by following two "best" values. The 
first one is the best solution (fitness) the particle has achieved so far. This value is called 
pbest. Another "best" value that is tracked by the particle swarm optimizer is the best value 
obtained so far by any particle within the neighborhood. This best value is a neighborhood 
or local best and called lbest. 

After finding the two best values, the particle updates its velocity and positions with 
following equations: 



a) 

(2) 



Wi is the inertia coefficient which slows velocity over time; v n is the current particle velocity; 
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p n is the current particle position in the search space; rl and r2 are random numbers 
between (0,1); cl and c2 are learning factors. The stop condition is usually the maximum 
number of allowed iterations for PSO to execute or the minimum error requirement. 

Because of the required search algorithm characteristics listed in section 1, we chose the PSO 
as the starting point for the search control algorithm. The PSO is computationally simple. It 
requires only four multiplies and four add/ subtracts to update the velocity and then one 
add to update the position. The PSO can also be a distributed algorithm. Each 
agent/ particle/ hot can update its own velocity and position. The only external information 
is the local best - the best value by any particle within the neighborhood. The calculation of 
the local best can be done with a simple comparative statement. Thus, each bot does not 
need to know the results from each member of the population as in many traditional 
schemes. (Though in some versions of the PSO, the lbest is replaced by gbest, the global 
best - the best value of any particle within the population.) 

The PSO also allows the contiguous movement of the bots. The updated position is relative 
to the current position so there are no jump changes in position or random movements. If 
there are constraints on the movement of the bot during each iteration, then limitations can 
be placed on the maximum and minimum velocity that is allowed for each particle/ bot. 

We propose embedding the PSO into a swarm of robots. In our approach, each bot is 
behaves as one particle in the PSO. Thus, each bot moves based on the PSO update 
equations (eqs. 1 and 2), makes a measurement and then broadcasts to the other bots in the 
swarm if it finds a new lbest measurement. We make some slight changes to the classic PSO 
algorithm to make it work with a robot swarm, so we call our algorithm the physically- 
embedded PSO (pePSO). 

Other authors have investigated using biological principles (Zarzhitsky & Spears, 2005; 
Valdastri et al., 2006; Schmickl & Crailsheim, 2006; Triannni et al., 2006) and PSO-type 
algorithms (Hayes et al., 2000; Hayes et al., 2002, Doctor et al., 2004; Pugh & Martinoli, 2006; 
Pugh & Martinoli, 2007; Jatmiko et al., 2006; Jatmiko et al., 2007) with multiple (simple) 
robots for search applications. Specifically, Hayes et al. report using autonomous mobile 
robots for beacon localization (Hayes et al., 2000) and plume tracking/ odor source 
localization (Hayes et al., 2002); they base their search techniques on biological principles 
(surge " upwind") but do not use the PSO algorithm directly. Doctor et al. (Doctor et al., 
2004) discuss using the PSO for multiple robot searches. Their focus is on optimizing the 
parameters of the search PSO and do not consider the scalability of the standard PSO to 
large numbers of robots. 

There are at least three other research groups that have investigated using mobile robots to 
do searches under the control of a PSO algorithm. Pugh et al. (Pugh & Martinoli, 2006; Pugh 
& Martinoli, 2007), explored using PSO on problems in noisy environments, focusing on 
unsupervised robotic learning. They used the PSO to evolve a robot controller that avoids 
obstacles and finds the source. They also investigated the possibility of using PSO without 
global position information. Jatmiko et al. (Jatmiko et al., 2006; Jatmiko et al., 2007) used 
mobile robots for plume detection and traversal. They utilized a modified form of the PSO 
to control the robots and consider how the robots respond to search space changes such as 
turbulence and wind changes. Akat and Gazi (Akat and Gazi, 2008) propose a version of 
the PSO for robot swarms that uses dynamic neighborhoods and asynchronous updates. 
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3. pePSO results 

3.1 Simulation Conditions 

We simulated the pePSO using three different test functions and five different target points 
with each function. The five different target points are shown in Figure 1. Only one target 
point was active during each simulation run; that is, there is only one target value in the 
search space at a time. 



Fig. 1. Map of simulation search space scaled to cover -1 to +1 and showing the 5 different 
target point locations. 


The three functions used in the simulation were the parabolic/ spherical function, Rastrigin 
function, and Rosenbrock function. The three test functions are given by: 


parabolic/ spherical: 

f(x,y) = (x-X lmget ) 2 +0- larger) 2 ( 3 ) 


Rastrigin: 

f(x,y) = (x- x target f - 10 cos(2^(x - x targel )) + 10) + (y - y tatgel f - 10 cos(2;r(y - y (argef )) + 10) 


Rosenbrock: 


f(x,y) = ( l-(x- 

•*7 arg et + 1)) 2 + 100«7 - y latget + 1) - (x 

X target + 1) 2 ) 2 (5) 

where (x,y) is the position of the bot/ particle and (x ta rget, y target) is the position of the target 
point. The spherical function was chosen because it approximates the expected dissipation 
pattern of chemicals, heat, etc that would be emitted by real-world targets. The Rastrigin 
and Rosenbrock functions were chosen because they are commonly used functions in PSO 
testing and they approximate a search space with obstacles and undulations. The 
Rosenbrock function used in our simulations is slightly different than the one reported in 
other simulations. We have shifted the minimum value of the function to the target location 
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instead of offset by (1,1) in x and y. This ensures that the target point is within the search 
space. 

Plots of the three test functions are shown in Figure 2. In each case, we are trying to find the 
location of the global minimum. The size of the search spaces were -100 to 100 for the 
parabolic function and -5.1 to 5.1 for the Rastrigin and Rosenbrock functions. (The 
dimensions correspond roughly to meters as we set the velocity on the bot based on meters 
per second.) Unlike most PSO simulations, the search space boundary for our simulations 
represents a "hard" border - the particles/ bots can not go outside the search space. We use 
a hard border because we want to approximate the conditions of actual robots searching in a 
room or some other confined space. The target point locations shown in Figure 1 were 
scaled for each of the search spaces. The (x,y) target locations are given in Table 1. 


Target 

point 

Parabolic 

Rastrigin/ 

Rosenbrock 

1 

(0,0) 

(0,0) 

2 

(94.7, 90.0) 

(4.86, 4.61) 

3 

(45.1, -98.5) 

(2.3, -5.02) 

4 

(-48.0, -52.7) 

(-2.46, -2.71) 

5 

(-79.6, 29.9) 

(-4.1, 1.54) 


Table 1. Target Locations for each of the three Test Functions 



Fig. 2. Plots of three test functions, (a) Parabolic function with contours. Min value (the 
target value) is at [-4.1 1.54]; (b) Rastrigin function with min (target) point at -2.46 -2.71; (b) 
Rosenbrock function with target point = [0 0]. Min value is at a saddle point. 
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The simulation incorporated the mobility limitations of our small robots. The max turning 
radius allowed was 36 degrees and the maximum velocity was 0.9 m/s. The bots were 
positioned randomly around the search space at t = 0. In our expected deployment scenario 
the bots will most likely be " dropped off" at a particular point (e.g., near a door or window) 
and not dispersed randomly about the search space. However, we can easily incorporate a 
dispersion algorithm (Hsiang et al. 2002; Morlok & Gini, 2004) to spread the bots 
throughout the search space before beginning the search phase using pePSO. 

For the simulation results, we used the following parameter values. Inertia coefficient, wi, 
was set to 0.9. The max velocity of the bots was set to 0.9 m/ s. Note that this is different 
than the old Vmax parameter in the original version of the PSO. We set an initial v for each 
bot/ particle to simulate the behavior of the physical robot. The cl and c2 coefficients were 
both set to 2.1. Since only three bots were used, the lbest topology was the same as gbest. 
That is, all of the bots communicated gbest and gbest location with all of the other bots. 

One major difference between our simulations and the PSO results reported in the literature 
is how we calculate a successful search, and, relatedly, the stop condition. A successful 
search occurred when a bot/ particle got within a tolerance of 0.2 in the x and y dimensions 
of the target location. Most PSO researchers track function value but we used location 
because of our search application. The simulation stopped when a bot got near the target or 
when 5400 function evaluations (about 15 minuntes) elapsed. Compared to other results, 
5400 is a relatively small number of function evaluations but we want the bots to find the 
target within a reasonable time period. 

3.2 Simulation results 

To evaluate the effectiveness of the pePSO, we made several simulation runs. An effective 
"hit" occurred when one of the bots found the global peak within 15 minutes of simulation 
time (5400 function iterations). We used the 5 different target points and did 200 test cases 
for each target point for a total of 1000 runs for each test function. In each test case, the 
initial start locations of the bots were different. We evaluated the overall effectiveness of the 
algorithm (i.e., how many times it found the "target") and compared the pePSO to the 
standard PSO. In tables 2, 3 and 4 we compare the pePSO to the standard PSO with the 
number of particles set to 3. We modified the standard PSO to calculate the percent found 
(number of hits) based on proximity to the target location rather than on the function value 
so that it was consistent with the pePSO. The tables report the percentage of targets found. 
The results are shown for each of the three test functions as well as the total for all five target 
points. 


For all three test functions, the pePSO performed much better than the PSO using same 
number of particles. 



Target point 

Target point 

Target point 

Target point 

Target 

Tot 


1 

2 

3 

4 

point 5 

al 

PSO 

100 

38.5 

51 

98 

90.5 

75.6 

pePSO 

99.5 

99.5 

99.5 

99.5 

100 

99.6 


Table 2. Parabolic Function, Results show the percent of targets found (out of 200 searches) 
for each of the five target points and also the overall results. 
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Target 
point 1 

Target 
point 2 

Target 
point 3 

Target 
point 4 

Target 
point 5 

Total 

PSO 

31 

33.5 

48 

26 

19 

31.5 

pePSO 

97.5 

96.5 

65.5 

99.5 

92 

90.2 


Table 3. Rastrigin Function. Results show the percent of targets foune (out of 200 searches) 
for each of the five target points and also the overall results. 



Target 
point 1 

Target 
point 2 

Target 
point 3 

Target 
point 4 

Target 
point 5 

Total 

PSO 

75 

53 

29.5 

51.5 

45 

50.8 

pePSO 

86.5 

92 

55.5 

88.5 

92 

82.9 


Table 4. Rosenbrock Function. Results show the percent of targets foune (out of 200 
searches) for each of the five target points and also the overall results. 


The results in Tables 2-4 show that none of the algorithms found the target location 100% of 
the time, even with a relatively easy problem like the parabolic function. There are a 
number of reasons for the missed detections. First, we declared a failed search after 
relatively few function evaluations (relative compared to other researchers' results); we 
limited both the standard PSO and the pePSO to 5400 function evaluations. Second, we are 
not using very many particles/ bots. Third, we have placed some of the target locations at 
positions near the edge of the search space which makes it somewhat more difficult for the 
PSO to find them. 

3.3 Hardware test conditions 

To investigate the effectiveness of the pePSO algorithm, we did several hardware 
experiments. In the experiments, a diffuse light source was placed near the ceiling of a dark 
room and pointed downward. Bots with light sensors were placed at various starting 
positions about 2 m apart to see (a) how often and (b) how quickly they could find the 
brightest spot of light in the room. 

The layout of the test area is shown in Figure 3. The hatched boxes in the middle represent the 
obstacles that we placed in the search space for the second half of the testing. The highest 
concentration of light is immediately to the left of the vertical obstacle. Since we are using a 
diffuse light source, the global best is actually a rectangle approximately .25 m by .3 m. 
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Fig. 3. Graphical representation of test area. Dimensions shown are in meters. Figure shows 
two obstacles and starred area is location of peak light intensity. 

To make the pePSO work in a robotic swarm, we had to make several adjustments. First, 
the bots determine their position by triangulating from three cricket motes set up as 
beacons. To correct for any missed packets, the bots are programmed to move to the next 
position and then wait for two consecutive "clean" measurements (distance measurements 
within 2 cm of each other) from all three beacons. This wait leads to relatively long search 
times. 

A second adjustment had to be made because of mechanical limitations in each bot. Since 
the bots steer with the front two wheels, the bots move in arcs rather than straight lines (see 
Figure 4). Each bot moves toward the desired position in the search space but upon arrival 
at the destination point, the orientation of the bot, as given by the direction of the wheels, is 
usually skewed relative to the movement. To compensate, we updated the bot's orientation 
angle at each iteration based on the bot's current position and its previous position. This 
eliminates the buildup of orientation errors that occur when a strict dead reckoning system 
for orientation is used. 



Fig. 4. Illustration of bot orientation mismatch. 
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A third adjustment was required because of the size of the bots. Unlike a simulation-only 
PSO, the hardware bots can get " stuck" at an obstacle or collide with another hot (particle). 
Once a hot got stuck or collided, we programmed the hot to back up and turn right. This 
allows the hot to move around long obstacles, such as a wall, even though it may require 
more than one cycle of backing up and turning to avoid. 

3.4 Hardware results 

During the hardware experiments, we programmed each hot in the swarm with identical 
programs. (The only difference is that each hot has a different identification number.) The 
PSO parameters used were cl = 2, c2 = 2, and Wi = 1.0. We tried using Wi = 0.8 but it slowed 
the bots down considerably and led to many failed searches. Each bot was programmed to 
move in the desired direction for approximately 0.5 sec. The bot would then make a 
measurement, determine its new position, calculate its desired movement direction based 
on the PSO update, orient its wheels to that direction, if possible, and then move for 0.5 sec. 
The search was ended when there had been 20 iterations of the algorithm with no new local 
best discovered. We made some test runs using a stop condition with 10 iterations with no 
new lbest but we determined that 10 was insufficient. 

Figure 5 shows the path traced out by 1 bot during a search sequence. The x and y axes are 
to scale and the axes are in cm. The asterisks (*) marks in the figures represent positions 
where a new local best was found. The rectangular box is the peak area of the search space. 
The bot starts in lower right corner. It moves up (north), finds new lbests and continues 
upward. When it starts to move away from the peak, it circles clockwise and then begins 
moving to the left (westward). When the light intensity measurements begin to taper off 
again, the bot circles counterclockwise back toward its previous best. The circle behavior is 
because the bot is limited in its turning radius - it can not make a sharp turn. Thus, it must 
move toward the global best in a roundabout fashion. Eventually, it settles on to the peak 
light value. 



Fig. 5. Path in search space with no obstacles for 1 bot. Asterisks show positions where lbests 
were recorded. 

The PSO algorithm was tested using swarm sizes of 1, 2 and 3 bots. (The 1-bot swarm was 
merely for baseline comparison.) Ten runs were made for each of the three cases. We 
tracked how many times the bots found the brightest spot and how long it took for the bots 
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to locate the peak. If twenty iterations of the algorithm passed with no new global best, then 
we declare the search is over. 

The quantitative results are shown in Tables 5 and 6. The results are for two different search 
spaces: one with no obstacles and one with obstacles. The table shows the number of bots in 
the swarm, the number of successful (completed) runs, the median search time for all 
searches, the average search time for just the successful searches, the standard deviation of 
the successful search times and the 90% confidence interval for the average time. The 
confidence interval is based on t-statistics using the successful searches as the degrees of 
freedom. 


Runs without obstacles. 


Test 

condition 

Complete 
runs (out of 
10) 

Median time 
(sec) 

Average 

(sec) 

Standard 
deviation (sec) 

Confidence 
Interval (90%) 

pePSO (1 
bot) 

6 

118.5 

180.2 

101 

±83.1 

pePSO (2 
bots) 

10 

177 

176.1 

68.9 

±39.9 

pePSO (3 
bots) 

10 

133.5 

109.6 

55.0 

±31.9 

Phototaxis 
(1 bot) 

9 

280 

273.8 

45.1 

±28.0 

3 bot (no 
commun) 

10 

351.5 

362.7 

98.7 

±57.2 


Table 5. Results from using PSO to program a suite of bots with sensor and position 
corrections. No obstacles in the search space 


Test 

condition 

Complete 
runs (out 
of 10) 

Median time 
(sec) 

Average 

(sec) 

Standard 

deviation 

(sec) 

Confidence 
Interval (90%) 

pePSO 

bot) 

(i 

8 

268 

250.2 

108.6 

±71.4 

pePSO 

bots) 

(2 

10 

178 

181.7 

76.1 

±44.1 

pePSO 

bots) 

(3 

10 

108.5 

125.8 

65.2 

±37.8 

3 bot 
commun) 

(no 

9 

290 

294.1 

111.3 

±69.0 


Table 6. Results from pePSO in a search space with obstacles 


Table 5 also shows the results from searches using a phototaxis approach and using three 
bots without any communication among the bots. These experiments provide a comparison 
for the effectiveness and search times for the pePSO. For the phototaxis experiments, we 
used a cylinder to make a directional light sensor on one of the bots. The bot then took light 
intensity measurements as it made a 360 degree loop. (The loop had a radius of 0.25 m.) The 
bot then moved in the direction of the greatest light reading. Our test results show that the 
multi-bot pePSO is faster and more effective than the phototaxis approach. 
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To compare the effectiveness of the multi-bot pePSO, we did a 3 hot search with no 
communication among the bots. Essentially, each hot moved toward bright regions of the 
room based on the PSO update equation without the lbest term. If a hot received a higher 
light reading than ever before, then it continued moving in that direction; thus, it is a 
pseudo-gradient method. The results in Table 6 show that the 3-bot pseudo-gradient 
method was effective (19 out of 20 successful runs) but it was much slower than the 3-bot 
pePSO. The 3-bot pePSO search is over 55% faster than the pseudo-gradient based 3-bot 
search both with obstacles and without obstacles. 

Increasing the number of bots does two things. It leads to more successful searches and 
reduces the time to find the peak/ best value. We see that even with the obstacles in the 
search space, the swarm is still able to find the " target" or peak light intensity every time for 
a multi-bot search. In general, the searches take longer when there are obstacles in the 
search space but the search is still successful. The times should be used for comparison 
purposes and not necessarily as an absolute reference. As mentioned earlier, sketchy 
communications with the beacons forced the bots to wait for two consecutive good data 
points at each iteration. This waiting time greatly increased the search time. 

To overcome the weak signal from one of the beacons, we modified the position algorithm 
to allow the bot to calculate its position with distance data from only two beacons. The bot 
initially tried to get two consecutive good measurements from all three beacons. If after 
seven seconds there was still not accurate distance information, then the bot would use the 
distance measurements from two beacons to calculate its (x,y) position. 

During our initial experiments, we noticed that the light sensors on the bots were 
mismatched. That is, different bots would read different values for the same light level. 
This imbalance led to some erroneous values for the search times. Specifically, if the bot 
with the lowest light readings found the peak value first, then other bots would circle 
toward that point. When another bot (with a light sensor that recorded higher light values) 
moved to the same location, it would record a higher light value and the algorithm would 
think a "new" target had been found. To correct for the different sensor readings, we used 
linear splines for each individual sensor to adjust and match the light sensor outputs of the 
three bots. 


4. Trophallaxis and robots 
4.1 Background 

In our trophallaxis-based algorithm, the bots (or, more generally, particles) do not actually 
exchange information but instead make sensor measurements when two or more 
bots/ particles are "in contact". The bots remain stationary for a certain time that is 
proportional to the measurement value. Bots thus cluster in areas of the search space that 
have high fitness/ measurement values. Each bot is independent (i.e., not reliant on neighbor 
bot measurements) but must have contact with the other bots to find the peaks. 

This new trophallaxis-based search algorithm has several advantages over other swarm- 
based search techniques. First, no bot-bot (particle-particle) communication is required. 
Thus, there is no concern with communication radius, protocol, or bandwidth. Unlike classic 
PSO-based techniques, the bots do not have to be arranged in topologies or communicate 
personal or global best information to any neighbors. The only communication that may be 
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required is signaling between bots to differentiate collisions between bots and collisions 
with obstacles. 

Second, the bots do not have to know their position. If position information is available, 
from beacons or some other source, then position information can be communicated at the 
end of the search. But during the search, the hot/ particle moves randomly except when it 
stops, takes a measurement, and waits. At the end of the search, the cluster locations can be 
determined from a remote camera, special-purpose robot, or human canvassing. 

Third, no on-board processing or memory is required - the bot does not even have to do the 
relatively simple PSO update equations. The bot/particle moves at random, takes a 
measurement and does a multiplication. It is so simple that a microcontroller may not be 
required, only some simple digital logic hardware. 

The Trophallactic Cluster Algorithm (TCA) has four basic steps: 

Step 1: Bots start randomly throughout the search space and then move at random through 
the search space. 

Step 2: If a bot intersects or collides with another bot, then it stops. 

Step 3: After stopping, the bot measures the " fitness" or function value at that point in 
space. It then waits at that point for a prescribed time based on the measurement. The 
higher the measurement value, then the longer the wait time. 

Step 4: When done, determine the locations of the clusters of bots. (We assume that this step 
is performed by an agent or agents that are separate from the swarm.) 

Step 1 is similar to the first step in the standard Particle Swarm Optimization (PSO) algorithm. 
For a software only optimization scheme, it is straightforward to randomly initialize the 
particles within the search boundaries. For a hardware scheme, a dispersion algorithm 
(Siebold & Hereford 2008; Spears et al., 2006) can be used to randomly place the bots. 

For random movement, we pick a direction and then have the bots move in a straight line in 
that direction until they encounter an obstacle, boundary or other bot. Thus, there is no path 
adjustment or Brownian motion type movement once the (random) initial direction is set. 
There is a maximum velocity at which the bots move throughout the search space. We 
experimented briefly with different maximum velocities to see its affect on overall results 
but we usually set it based on the expected maximum velocity of our hardware bots. 

In step 2, we detect collision by determining whether bots are within a certain distance of 
each other. In software, this is done after each time step. In hardware, it can be done using 
infrared sensors on each bot. 

In a hardware implementation of the TCA, we would need to distinguish between collisions 
with obstacles and collisions with other bots. Obstacles and walls just cause the bot/ particle 
to reorient and move in a new direction. They do not lead to a stop/ measure/ wait 
sequence. So once a collision is detected, the bot would have to determine if the collision is 
with another bot or with an obstacle. One way to do this is to have the bots signal with 
LEDS (ala trophallaxis where only neighbors next to each other can exchange info). In this 
way, each bot will know it has encountered another bot. 

Once a bot is stopped (as a result of a collision with another bot), then it measures the value 
of the function at that location. (In hardware, the bot would take a sensor reading.) Since the 
bots are nearly co-located, the function value will be nearly the same for both bots. The wait 
time is a multiple of the function value so the bot(s) will wait longer in areas of high fitness 
relative to areas with low function values. Other bots may also "collide" with the stopped 
and waiting bots but that will not reset the wait times for the stopped bots. For the results in 
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this paper the wait time was exponentially related to the measurement value; we 
experimented with linear wait times as well. 

We did step 4 (determine the clusters) when the search is "done". In general, the bots begin 
to collide/ stop/ wait at the beginning. Thus, the bots tend to cluster soon after the search 
begins so the search can be stopped at any time to observe the location(s) of the clusters. In 
2D, the clusters tend to become more pronounced as the iterations increase so waiting a 
longer time can make the position(s) of the peak(s) more obvious. 

4.2 Related work 

The TCA is based on the work of Thomas Schmickl and Karl Crailsheim (Schmickl & 
Crailsheim 2006; Schmickl & Crailsheim 2008) who developed the concept based on the 
trophallactic behavior of honey bees. Schmickl and Crailsheim use the trophallactic concept 
to have a swarm of bots move (simulated) dirt from a source point to a dump point. The 
bots can upload "nectar" from the source point, where the amount of nectar for each bot is 
stored in an internal variable. As the robots move, the amount of stored nectar decreases, so 
the higher the nectar level, then the closer to the source. Each robot also queries the nectar 
level of the robots in the local neighborhood and can use this information to navigate uphill 
in the gradient. There is a also a dump area where the loaded robots aggregate and drop the 
"dirt" particles. The swarm had to navigate between the source and the dump and achieved 
this by establishing two distinct gradients in parallel. 

Their preliminary results showed a problem where the bots tended to aggregate near the 
dump and the source. When that happened, the bridge or gradient information between the 
source and the dump was lost. To prevent the aggregation, they prevented a percentage of 
their robots from moving uphill and just performed a random walk with obstacle avoidance. 
Even though the work of Schmickl and Crailsheim is significant, they show no published 
results where they apply the trophallactic concept to strictly search/ optimization type 
problems. Nor do they show results when there is more than one peak (or source point) in 
the search space. They also require bot-bot communications to form the pseudo-gradient 
that the loaded (or empty) bots follow, while our TCA approach does not require adjacent 
particles/ bots to exchange nectar levels (or any other measured values). 

In (Ngo & Schioler, 2008), Ngo and Schioler model a group of autonomous mobile robots 
with the possibility of self-refueling. Inspired by the natural concept of trophallaxis, they 
investigate a system of multiple robots that is capable of energy sharing to sustain robot life. 
Energy (via the exchange of rechargeable batteries) is transferred by direct contact with 
other robots and fixed docking stations. 

In this research, we are applying the concept of trophallaxis to solve a completely different 
type of problem than Ngo and Schioler, though some of their results may be applicable if we 
expand our research to include energy use by the robots. 


5. Trophallaxis search results 
5.1 Test conditions 

We tested the TCA algorithm on three functions: two ID and one 2D. The ID functions are 
denoted F3 and F4 and were used by Parrott and Li (Parrott and Li, 2006) for testing PSO- 
based algorithms that find multiple peaks in the search space. The equations for F3 and F4 
are given by 
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F3(x) = sin 6 (5;r(x 3/4 - 0.05)) (6) 

r _ A AO 

F4(x) = exp(-21og(2)(— — — ) 2 )sin 6 (5^(x 3/4 - 0.05)) (7) 

.o j4 

Plots for the function F3 and F4 are shown in figure 6. Each ID test function is defined over 
the scale of 0 < x < 1. F3 has five equal-height peaks (global optima) that are unevenly 
spaced in the search space. F4 also has five unevenly spaced peaks but only one is a global 
optimum while the other four peaks are local optima. Our goal is to find all five peaks; that 
is, the global optimum plus the local optima. The peak locations are given in Table 7. 




X 


Fig. 6. ID test function with equal peaks, F3, and with unequal peaks, F4 


Peak# 

1 

2 

3 

4 

5 

X locations 

.0797 

.2465 

.4505 

.6815 

.9340 


Table 7. Peak locations for test functions F3 and F4 

The 2D function is a slight variation of the standard Rastrigin function. The equation for the 
Rastrigin is given in equation 4 with the range on x and y is -5.12 to 5.12. The Rastrigin is 
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highly multimodal (see figure 2) and has one global minimum. For the TCA simulations, 
we modified the Rastrigin so that it has a peak of 1 at (.35, .35) instead of a minimum at the 
origin. We also scaled the function slightly so that there are nine peaks (1 global, 8 local) 
within the [-5.12, 5.12] range. The modified Rastrigin function is shown in figure 7. 



Fig. 7. Plot of modified Rastrigin function scaled so that it has 9 peaks in [-5.12,5.12] range 
and peak value is equal to 1. 

We evaluated the effectiveness of the TCA algorithm using two different metrics. The first 
metric is the total percentage of peaks found (found rate). Since each function has multiple 
peaks (both global and local) we totaled the number of actual peaks that the swarm found 
divided by the total number of peaks (see equation 8). Note that "peaks found" refers to 
only those hot clusters that are within ± .04 (ID) or a radius of .4 (2D) of the actual peak 
location. 

Found rate = (peaks found)/ (total number of peaks in search space) (8) 

The second metric is related to the success rate from Parrott and Li: 

Success rate = (Runs where more than half of total peaks found)/ (total number of runs) (9) 
The success rate gives an indication of how many of the runs were successful. We designate a 
successful run as one where a majority of the peaks (more than half) in the search space are 
found. This definition is based on the robot search idea where the goal is for the robots to 
quickly cluster near the target points. We note that this is a slightly different definition for 
success rate than Parrott and Li; their success rate is based on the closest particle and finding 
all of the peaks in the search space. 

The locations of the bot clusters were determined using the K-means clustering algorithm. 
The K-means algorithm minimizes the sum, over all clusters, of the point-to-cluster centroid 
distances. We then compared the cluster centroid to the actual peak location to determine if 
the peak was found or not. For the ID functions, we used a tolerance of ± 0.04 between the 
cluster centroid and the actual peak and for the 2D functions we used a radius of 0.4. 

We considered two clustering approaches; the first one uses all of the bots when 
determining the cluster centroids. The second approach uses only the final position of the 
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bots that are stopped (that is, in a collision) when determining clusters. We refer to the 
second approach as " cluster reduction", since it reduces the number of bots that are 
considered. 

5.2 Trophallactic Cluster Algorithm ID results 

Qualitative results from computer simulations of the TCA for the two ID functions F3 and 
F4 are shown in Figure 8. The top plot in the figure shows the original function; the middle 
plot shows the final bot positions (after 400 iterations) with each bot position represented by 
a star (*). The bottom plot is a normalized histogram; the histogram is made by tracking the 
position of each bot after each time interval. The figure reveals that the bots do cluster 
around the peaks in the function and thus give evidence that the TCA will reliably find 
multiple peaks. 

The histogram plots reveal some interesting information. For both F3 and F4, there are 
significant peaks in the histogram at the same locations as the function peaks, providing 
evidence that the bots spend a majority of time near the peaks and it is not just at the end of 
the simulation that the bots cluster. Also, for F3 the histogram peaks have approximately the 
same amplitude (peak amplitudes range from 0.7 to 1.0). For F4, however, the histogram 
peaks diminish in amplitude in almost direct proportion to the function peaks (peak 
amplitudes diminish from 1.0 down to 0.25). This implies that the bots are spending more 
time in the vicinity of the larger peaks. The bots are thus "attracted" to stronger signals, as 
expected. 
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Fig. 8. Qualitative results for function F3 (left) and function F4 (right). Final bot locations are 
shown in the middle plot and histograms of bot positions are shown in the bottom plot 


We performed computer simulations to tailor three of the parameters of TCA algorithm for 
ID functions. The three parameters were tmax, the maximum number of iterations for the 
simulation, nbots, the number of bots to use, and waitfactor. The waitfactor sets how long each 
bot waits based on the measured value after a collision. We tried linear wait functions (wait 
time increases linearly with measurement value) but had more success with exponential 
wait functions give by 
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Wait time = waitfactor * (e( measurement )-l) 


(ii) 


For the parameter selection, we varied one parameter at a time and repeated the simulation 
100 times. Plots of the average found rate for parameters nbots and waitfactor with no cluster 
reduction are shown in Figure 9. The first plot shows the found rate as nbots is varied from 
10 to 200 with tmax set to 500 and waitfactor set to 5. The second plot shows the found rate as 
waitfactor is varied from 1 to 10 with tmax = 500 and nbots = 80. Similar tests were done with 
cluster reduction. 



0 20 40 


£ 0.75 


60 80 100 120 140 160 180 200 

Number of Bots 


i r - 


Wait Factor 

Fig. 9. Results showing found rate vs nbots (left figure) and waitfactor (right figure) for F3 
function 


When the peaks were found with no cluster reduction, the found rate versus the parameter 
value curve resembled (l-e- x ) shape and asymptotically approached a found rate of about 
78%. Thus, there was not one precise parameter value but a range of parameter values that 
led to the best performance: nbots greater than 80 and waitfactor greater than 3. The tmax 
curve was flat - it appears that the bots quickly cluster near the peaks and there is little 
change in performance as the tmax is increased. A summary for the parameter selection 
process is shown in Table 7. 
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Parameter 

w/out cluster reduction 

w/ cluster reduction 

nbots 

>80 

<20 

tmax 

>300 

>100 

waitfactor 

>3 

>4 


Table 7. Best parameter ranges for TCA for ID functions 


When the hot cluster centroids were found with cluster reduction, the response curves for 
tmax (flat) and waitf actor (exponential) were similar in shape as without cluster reduction, 
though the curve asymptotically approached a found rate of 86%. The response curve for 
nbots was different, however. The found rate went up as nbots decreased so fewer bots was 
better than more bots, assuming that there at least 5 stopped bots in the search space. It 
appears that for a small number of bots, that there was a smaller percentage of bots in the 
clusters (say 13 out of 20 instead of 85 out of 90) and those off-clusters bots moved the 
cluster centers away from the peaks and led to missed detections. 

For the final results we used the parameter values tmax = 500, waitfactor = 4, and nbots = 60. 
We used the same parameter values for all test cases: two different functions (F3 and F4) 
and with and without cluster reduction. There was a slight dilemma on the choice for nbots 
since more bots did better without cluster reduction and fewer bots did better without 
cluster reduction so we compromised on 60. We ran the ID simulations 500 times and the 
results are shown in Table 8. 



F3 

F4 


Avg # 

peaks 

found 

Std 

dev 

Found 

rate 

Success 

rate 

Avg # 

peaks 

found 

Std 

dev 

Found 

rate 

Succes 

s rate 

w/out 

cluster 

reduction 

4.052 

.9311 

81.04 

% 

97.2% 

4.016 

.9195 

80.32% 

97.6% 

w/ cluster 
reduction 

4.130 

.7761 

82.6% 

98.2% 

4.098 

.7782 

81.96% 

99.0% 


Table 8. Final results showing average number of peaks found (out of 5 peaks), found rate 
and success rate for 500 iterations 


The ID results show that the TCA was very effective at finding a majority of the peaks in the 
2 different functions. The success rate was above 97% and the found rate was above 80%. 
These are good results for an algorithm where the individual particles/ bots do not have 
position information and no bot-bot communication is required. 

The results shown in Table 8 are very consistent. The TCA algorithm finds 4 out of the 5 
peaks and there is little difference in the results between F3 (peaks of equal height) and F4 
(peaks have different heights). There is a slight improvement with cluster reduction, that is, 
when only the stopped bots are used to determine the peak locations. 

5.3 Trophallactic Cluster Algorithm 2D results 

The results of a typical two dimensional search using the TCA are shown in Figure 10. The 
first figure shows a plot of the Rastrigin function with the final bot positions superimposed 
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on top of it. The second figure shows the final hot positions and the centroids of nine 
clusters found by the K-means clustering algorithm (black stars). Note that six of the cluster 
centroids are close to actual peaks in the Rastrigin function, but only one of the centroids 
was within the required tolerance and was thus declared a peak (red diamond). 

The initial 2D results, like those shown in Figure 10, illustrate the usefulness of cluster 
reduction. Figure 11 shows the same final bot positions as in Figure 10, except only the 
clusters with three of more bots are kept. That is, small clusters of two bots and any bots not 
in a cluster are eliminated. The K-means clustering is performed with this smaller set of bots 
and the cluster centroids compared to the peak locations. 

After cluster reduction, there are four cluster centroids that are within the tolerance radius 
of the peak instead on only one centroid. Thus, ignoring the still-moving bots after the 
conclusion of the search clarifies the definition of the clusters of bots. This in turn leads to 
more accurate identification of the peaks in the function. 

As with the ID test functions, computer simulations were conducted to refine the three 
parameters of tmax, nbots, and waitfactor for the 2D case. The results from these simulations 
for nbots and tmax are shown in Figure 12. Each graph shows the found rate as the 
parameter was varied; they are the result of averaging 100 simulations for each set of 
parameters. For found rate vs nbots graph, tmax was set to 1600 and waitfactor was set to 4. 
For the found rate vs tmax graph, nbots= 300 and waitfactor = 4. 
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Fig. 10. Typical TCA search results for the Rastrigin 2D function. Left figure: Rastrigin 
function showing final hot position. Right figure: final hot position with cluster centroids - 
red diamond denotes found peak and black star denotes cluster centroid. 
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Fig. 11. Analysis of typical TCA search with cluster reduction. Red diamond denotes found 
peak. Black star denotes inaccurate peak. 
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The results and interpretation of these two dimensional results are similar to the one 
dimensional case. The results roughly follow a (l~e x ) form. Therefore, the appropriate 
parameter values are again ranges rather than precise values. The values are given in Table 9. 
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Fig. 12. Results showing found rate vs nbots (left figure) and tmax (right figure) for Rastrigin 
function 



Parameter 

w/out cluster reduction 

w/ cluster reduction 

nbots 

>500 

<300 

tmax 

>1600 

>1700 

waitf actor 

>4 

>4 


Table 9. Best parameter ranges for 2D Rastrigin function 


The final two dimensional results were obtained using parameter values tmax = 1600, nbots 
= 600 , audio aitf actor = 4. The same parameters were used both with and without cluster 
reduction. We averaged the results from 500 simulations and the results are shown in Table 
10. 
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Avg # peaks found 

Std deviation 

Found rate 

Success rate 

w/out cluster reduction 

3.7360 

1.4375 

41.5% 

29.2% 

w/ cluster reduction 

3.3820 

1.4888 

37.6% 

21.8% 


Table 10. Final results showing average number of peaks found (out of 9 peaks), found rate 
and success rate for 500 iterations for the Rastrigin 2D function 


The results from the 2D Rastrigin function are not as good as the results from the ID 
functions. The lower found rate is due primarily to the fact that the Rastrigin function is a 
hard function - the peaks do not stand out as prominently as the F3 or even the F4 peaks. In 
addition, the 2D search space is much larger; for the Rastrigin function, we used a scale of - 
5.1 to +5.12 for both x and y, while the ID functions are only defined between 0 < x <1. We 
increased the tolerance for the 2D results to 0.4 and it appeared that many cluster centroids 
were close to the actual peaks but, unfortunately, not within the tolerance radius. 

6. Conclusions 

We developed and tested two biologically inspired search strategies for robot swarms. The 
first search technique, which we call the physically embedded Particle Swarm Optimization 
(pePSO) algorithm, is based on bird flocking and the PSO. The pePSO is able to find single 
peaks even in a complex search space such as the Rastrigin function and the Rosenbrock 
function. We were also the first research team to show that the pePSO could be 
implemented in an actual suite of robots. 

Our experiments with the pePSO led to the development of a robot swarm search strategy 
that did not require each bot to know its physical location. We based the second search 
strategy on the biological principle of trophallaxis and called the algorithm Trophallactic 
Cluster Algorithm (TCA). We have simulated the TCA and gotten good results with multi- 
peak ID functions but only fair results with multi-peak 2D functions. The next step to 
improve TCA performance is to evaluate the clustering algorithm. It appears that many 
times there is a cluster of bots near a peak but the clustering algorithm does not place the 
cluster centroid within the tolerance range of the actual peak. A realistic extension is to find 
the cluster locations via the K-means algorithm and then see if the actual peak falls within 
the bounds of the entire cluster. 
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1. Introduction 

Nowadays the level of automation in manufacturing industries has been increased 
dramatically. Some examples of these automation progresses are in cellular manufacturing 
and robotic cells. A growing body of evidence suggests that, in a wide variety of industrial 
settings, material handling within a cell can be accomplished very efficiently by employing 
robots (see (Asfahl, 1992)). Among the interrelated issues to be considered in using robotic 
cells are their designs, the scheduling of robot moves, and the sequencing of parts to be 
produced. 

Robotic cell problem in which robot is used as material handling system received 
considerable attentions. Sethi et al. (1992) proved that in buffer-less single-gripper two- 
machine robotic cells producing single part-type and having identical robot travel times 
between adjacent machines and identical load/unload times, a 1-unit cycle provides the 
minimum per unit cycle time in the class of all solutions, cyclic or otherwise. For three 
machine case, Crama and van de Klundert (1999), and Brauner and Finke (1999) shown that 
the best 1-unit cycle is optimal solution for the class of all cyclic solutions. Hall et al. (1997; 
1998) considered the computational complexity of the multiple-type parts three-machine 
robotic cell problem under various robot movement policies. This problem is studied for no- 
wait robotic cells too. For example Agnetis (2000) found an optimal part schedule for no- 
wait robotic cells with three and two machines. Agnetis and pacciarelli (2000) have studied 
partscheduling problem for no-wait robotic cells, and found the complexity of the problem. 
Crama et al. (2000) studied flow-shop scheduling problems, models for such problems, and 
complexity of theses problems. Dawande et al. (2005) reviewed the recent developments in 
robotic cells and, provided a classification scheme for robotic cells scheduling problem. 
Some other special cases have been studied such as: Drobouchevitch et al. (2006) provided a 
model for cyclic production in a dual-gripper robotic cell. Gultekin et al. (2006) studied 
robotic cell scheduling problem with tooling constraints for a two-machine robotic cell 
where some operations can only be processed on the first machine and some others can only 
be processed on the second machine and the remaining can be processed on both machines. 
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Gultekin et al. (2007) considered a flexible manufacturing robotic cell with identical parts in 
which machines are able to do different operations and the operation time is not system 
parameter and is variable. They proposed a lower bound for 1-unit cycles and 2-unit cycles. 
Sriskandarajah et al. (1998) classified the part sequence problems associated with different 
robot movement policies, in this chapter a robot movement policy is considered, which its 
part scheduling problem is NP-Hard, and Baghchi et al. (2006) proposed to solve this 
problem, by a heuristic or meta-heuristic. In this chapter a meta-heuristic method based on 
particle swarm optimization is applied to solve the problem. 

In this chapter an m-machine flexible cyclic cell is considered. All parts in an MPS (A 
minimal part set) visit each machine in the same order, the production environment is 
cyclic, and parts are produced at the same order repeatedly. 

In this chapter, we consider multiple-type parts three-machine robotic cells which have 
operational flexibility in which the operations can be performed in any order; moreover 
each machine can be configured to perform any operation. To explain the problem, consider 
a machining centre where three machine tools are located and a robot is used to feed the 

machines namely M x JVf 2 J\4 2 (see figure 1). All parts are brought to and removed from the 

robotic cell by Automated Storage & Retrieval System (AS/RS). The pallets and feeders of 
the AS/RS system allow hundreds of parts to be loaded into the cell without human 



Robot 


AS/RS 

Fig. 1. Robotic work cell layout with three machines 

The aim of this chapter is to find a schedule for the robot movement and the sequence of 
parts to maximize throughput (i.e., to minimize cycle time), as it is showed that this problem 
is NP-Complete in general (see Hall et al. (1997)). Hence, this chapter proposes a novel 
hybrid particle swarm optimization (HPSO) algorithm to tackle the problem. To validate the 
developed model and solution algorithm, various test problems with different sizes is 
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randomly generated and the performance of the HPSO is compared with three benchmark 
metaheuristics: Genetic Algorithm, PSO-I (basic Particle Swarm Optimization algorithm), 
and PSO-II (constriction Particle Swarm Optimization algorithm). The rest of this chapter is 
organized as follows: The problem definition and required notations are presented in 
Section 2, Section 3 presents the developed mathematical model, and in Section 4, the 
proposed hybrid particle swarm optimization algorithm is described. The computational 
results are reported in Section 5, and the conclusions are presented in Section 6. 


2. Problem definition 

The robotic cell problem is a special case of the cyclic blocking flow-shop, where the jobs 
might block either the machine or the robot. In a cyclic schedule the same sequences repeat 
over and over and the state of the cell at the beginning of each cycle is the similar to the next 
cycle. It is assumed that the discipline for the movements of parts is an ordinary flow-shop 

discipline. That is a part meets machines M x J\P 2 JVt 3 consequently. 


2.1 Notations 

The following notation is used to describe the robotic cell problem: 

TTl : The number of machines 

I/O : The automated input-output system for the cell 

PT f : The part-type i to be produced 

r. : The minimal ratio of part i to be produced 

MPS : The number of part set consisting r. parts of type PP 

Tl : the total number of parts to be produced in the MPS ( Tl — V x + V 2 + . . . + r k ) 

CL i : The processing time of part i on M x 

b { : The processing time of part i on M 2 

C { : The processing time of part i on M 3 

: Robot travelling time between two successive machines (I/O is assumed as 
^ machine M 0 ) 

S : The load/ unload time of part i 

wj : The robot waiting time on Mj to unload part i 

$ k : The robot movement policy S under category k 
T k : The cycle time under S k 

In this study the standard classification scheme for scheduling problems: \f/ x \ lf/ 2 \ y/ 2 is 
used where l// l indicates the scheduling environment, lf/ 2 describes the job characteristics 
and Iff 3 defines the objective function (Dawande et al., 2005). For example 
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FRC 3 | k > 2, S l | C t denotes the minimization of cycle time for multi-type part problem 
in a three flow-shop robotic cell, restricted to robot move cycle S l . 

2.2 Three machine robotic flow shop cell FRC 3 \ K > 2 \ C t 

In the three machine robotic flow shop cell, there are six different potentially optimal 
policies for robot to move the parts between the machines (Bagchi et al., 2006). Sethi et al. 
(1992) showed that any potentially optimal one-unit robot move cycle in a m machine 
robotic cell can be described by exactly m+1 following basic activities: 

: Load a part on 1 — 1, 2 , . . TYl 

M+ : Unload a finished part from i = 1, 2 , . . TYl 

In other words, a cycle can be uniquely described by a permutation of the m+1 activity. The 
following are the available robot move cycles for m=3 flow-shop robotic cell (Sethi et al., 
1992): 

5' : {M3 + , My , M ~ , M ~ , Ml } 

V : [m; , My , m ~ , m~ , m: } 

S 3 : [m; , M ~ , My , M~ , m; } 
s 4 -.{m;,m~, m~ , M- , m: } 

S 5 :{m;,M-, M y , M~ , m : } 

S 6 : [m; , M~ , M~ , My , m; } 

In this chapter we consider a three machine robotic cell problem under the S 6 policy 

(Figure 2). The problem of finding the best part sequence using the robot move cycle S 6 is 
NP-complete (Hall et al., 1998). 
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Fig. 2. The robot movement under S 6 
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Lemma 1 . The cycle times of one unit for the policy S are given by: 


,6 



Pickup part p i+2 from I/O (s) move it to Mj (S) load p i+2 onto M 1 (s) go to M 3 (2 §) if 
necessary wait at M 3 (wh) , unload p t from M 3 (s) move it to I/0(£) drop p t at I/0(^) 
go to M 2 (2 5) if necessary wait at M 2 (w 1+1 2 ) , unload p i+l from M 2 (s) move it to M 3 (S) , 

load P i+l onto M 3 (s) go to Mj (2 S) if necessary wait at M^w 14 ^), unload P i+2 from 
Mj {s) move it to M 2 (S) load P i+2 onto M 2 (s) go to I/O (IS) then start a new cycle by 
picking up the part P i+3 . 

The cycle time by considering waiting times is as follow: 



= \2S + 8s + w[ +2 + Wp + W 3 


I ,<j(i)cr(i+\)<j(i+2) 



3. Developing mathematical model 

In this section we develop a systematic method to produce necessary mathematical 
programming formulation for robotic cells. Therefore first we model single-part type 
problem through Petri nets, and then extend the model to multiple-part type problem. 

A Petri-net is a four-tuple PN(P , T , A ? W) , where P = {p x , p 2 , p n } is a finite set of 
places, T = {t l ,t 2 ,...,t m } is a finite set of transitions, A (Z (Px r)u(rxP) is a finite 
set of arcs, and W \ A {1, 2, 3, ... } is a weight function. 

Every place has an initial marking M 0 P — > {0,1, 2,...} . If we assign time to the 
transitions we call it as Timed Petri net. 

The behaviour of many systems can be described by system states and their changes, to 
simulate the dynamic behaviour of system; marking in a Petri-net is changed according to 
the following transition (firing) rule: 1) A transition is said to be enabled if each input place 
p of f is marked at least with t ) tokens, where w(j^, t) is weight of the arc from p to 

t. 2) An enabled transition may or may not be fired (depending on whether or not the event 
takes place). A firing of an enabled transition t removes w(p,t) tokens from each input 
place p of t and adds w(p,t) tokens to each output place p of t , where w(p,t) is the 
weight of the arc from t to p. 
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By considering a single-part type system, the robot arm at steady state is located at machine 
M 2 , therefore by coming back to this node we have a complete cycle for the robot arm. 

The related Petri net for robot movements is shown in Figure 3 and the descriptions of the 
nodes for this graph with respective execution times would be as follows: 



Fig. 3. Petri net for S 6 policy 


R x : go to M 3 ((5) ; R 2 : load M 3 (s) ; 

R 4 : unload M x (g) ; R 5 : go to M 2 (25) ; 

R 7 : go to input, pickup a new part, go to M x (s + 35 ) ; 

R 9 : go to M 3 (25) ; R xo : unload M 3 (s) ; 

R xx : go to output, drop the part, go to M 3 (6* + 3(5) ; 

7?/^ : wait at M . (w l j ) S- : starting time of 7?. ; 
a : Mj is ready to be unloaded; 

P : M 2 is ready to be unloaded; 

Y : M 3 is ready to be unloaded; 

By considering a multiple-part type system, at machine M x , when we want to load a part 
on the machine we have to decide which part should be chosen such that the cycle time is 


7? 3 : go to M x (25) ; 

7? 6 : load M 2 (s) ; 

7? 8 : load M x (<?); 

R n : unload M 2 ( 6 *) ; 
sp : starting time of RP 
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minimized. The same thing also can be achieved for M 2 and M 3 . Based on the choosing 
gate definition we simply have three choosing gates as CC , (5 , and y . Thus we can write 
the following formulation using 0-1 integer variables xly , x2 { . , and x3 fj . as: 

n 

«1 ' S A,X ~ \n + C t ^ E Xl in ( U i ) + e 

i = 1 
n 

a j :s 4j + i- s sj^'Z xl y( a i) + £ j = 2,-,n. 

i= 1 
n 

P] :S i2,J- S 6,J ^Ya X2 ip b P + S j = h-,n. 

i= 1 
n 

Yj - s w ,j- s 2j^Yj x2> iP c P + £ j = 

i= 1 

Definition. A marked graph is a Petri-net such that every place has only one input and only 
one output. 

Theorem 1. For a marked graph which every place has i m tokens (see figure 4), the 
following relation S B > S A + 7! ijC t , where S A , S B are starting times of transitions A and B 

respectively, and C t is cycle time, is true. 

A p B 

□ — ►© — 

Fig. 4. The marked graph in theorem 1 
Proof: see ref. (Maggot, 1984). 

In addition the following feasibility constraints assign unique positioning for every job: 

n 

i = i 
n 

Yj x h = l i = 

y=i 

To keep the sequence of the parts between the machines in a right order, we have to add the 
following constraints: 

xl,j =x2. J+1 i = 1 = 1 ,—,n 

xl Uj =x3 iJ+1 i = 1 = 1, •••,«. 

Where, we assume that xl i n+1 = xl i X because of the cyclic repetition of parts. 

Thus the complete model for the three machine robotic cell with multiple-part would be as 
follows: 
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min Cf 


Subject to: 


Pl, 1 ^2,1 -*12,1 + C t + 8 

j = V 

• ,n 

a) 

Pl,j ■ S 2j ~ S 12j = £ + 8 

7 = 1," 

■ ,n 

(2) 

P3,j :S 4j ~ S 2 j~ W lj =£ + 2S 

7 = 1," 

■ ,n 

(3) 

P5, j ■ s 6 j ~ s 4j = £ 4- S 

7=1," 

-,n 

(4) 

Pl ,j • s 8 j ~ s 6 j = 2 ^ + 3S 

7 = 1," 

‘,n 

(5) 

P9,j '■ s 10 j ~ s Sj ~ w 3j = e + 25 

7 = 1," 

■ ,n 

(6) 

P\ 1 J ■H7.j-s\0j-W2j= 2 £ + 3^ 

7 = 1," 

-,n 

(7) 

Yl 

• *^4,1 — ^7,1 ~ x ^in ( a i *P aj ) — 8 

i=l 



(8) 

n 

<x ] • s 4,j ~ s l,j ~ ^ x ^ij^ a i VcLf ) — 8 
i= 1 

7 = 2,- 

• ,n 

(9) 

n 

P] : S \2J ~ S 6J ~ ^ x2 ij( b i + <Pb t ) ^ £ 

i=l 

7 = 1," 

■ ,n 

(10) 

n 

Y]' S WJ ~ s 2 J-^ x3 ij( c i + ^P Ci ) * 8 

i=l 

x Kj-i =x2 U 

7 = 1," 

■ ,n 

(11) 

U 7=1,- 

•• ,n 

(12) 

x2 i,j- 1 = x2> i,j 

1,7 = 1,- 

•• ,n 

(13) 

II 

Hi 

7 = 1," 

■ ,n 

(14) 

n 

y 1 / 1 /'/ = 1 

1 = 1," 

-,n 

(15) 

7=1 




S i j > 0, W k j > 0, xl, x2, x3 e {0,1} 





4. The proposed hybrid particle swarm optimization (HPSO) algorithm 

The particle swarm optimization (PSO) is a population based stochastic optimization 
technique that was developed by Kennedy and Eberhart in 1995 (Hu et al., 2004). The PSO 
inspired by social behavior of bird flocking or fish schooling. In PSO, each solution is a bird 
in the flock and is referred to as a particle. A particle is analogous to a chromosome in GAs 
(Kennedy and Eberhart, 1995). All particles have fitness values which are evaluated by the 
fitness function to be optimized, and have velocities which direct the flying of the particles. 
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The particles fly through the problem space by following the particles with the best 
solutions so far (Shi and Eberhart, 1998). 

The general scheme of the proposed HPSO is presented in Figure 5. 



Fig 5. The schematic structure of the proposed HPSO 


In this chapter, we extend the discrete PSO of Liao et al. (2007) to solve the robotic cell 
problem. In the proposed HPSO the velocity of each particle is calculated according to 
equation (16). 

V id = w V id + q x rand () x ( pBest id - x id ) + c 2 x Rand () x ( lBest id - x id ) - ( 16 ) 
it /(Number of Iterations )) x(ax Frequency Matrix - b ) 


Where C } and C 2 are the learning factors that control the influence of pBest and IBest. W is 
the inertia weight which controls the exploration and exploitation abilities of algorithm. 
rand () and Rand () are two independently generated random numbers, t is the current 

iteration and a and b are two parameters that adjust the influence of the Frequency 
Matrix on velocity value. pBest is the best position which each particle has found since the 
first step and it represents the experiential knowledge of a particle. After the cloning 
procedure (the detailed of cloning procedure will be described in the next section), a 
neighborhood for each particle is achieved. The best particle in this neighborhood is selected 
as IBest. 
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As Liao et al. (2007), the velocity values transfers from real numbers to the probability of 
changes by using the equation (17): 


(17) 



where s(V id ) stands for the probability of X id taking the value 1. In the proposed 

algorithm, the new position (sequence) of each particle is constructed based on its 
probability of changes that calculated by equation (17). Precisely, for calculating the new 
position of each particle, the algorithm starts with a null sequence and places an 
unscheduled job; in position k (k= 1, 2, n) according to the probability that determined 
by equation (18): 



(18) 


where F is the set of the first f unscheduled jobs as present in the best particle (solution) 
obtained till current iteration. To achieve a complete sequence, the jobs are added one after 
another to the partial sequence. 

The proposed HPSO terminates after a given number of iterations and the best sequence is 
reported as the final solution for the problem. 

4.1 Cloning 

For avoiding local optimal solutions we implement cloning procedure which in summary 
can be described as follows: 

1. M copies of the solution are generated so that there are (M+l) identical solutions 
available. 

2. Each of the M copies are subjected to the swapping mutation. 

3. In each clone only the original solution participates in HPSO evolution procedure 
whereas the other copies of the solution would be discarded. 

4. The above procedure is repeated for all of the solutions in the swarm. 

4.2 Fitness evaluation 

As any metaheuristic algorithm, the HPSO uses a fitness function to quantify the optimality 
of a particle (sequence). The cycle times of one unit for the policy S 6 are given by: 



= \25 + ^>e + max{0,a CT(i+2) -8£-4£,b CT(i+1) -8c>-4£,c CT(i) -8 5 As} 


I,cr(i)cr(i+l)cr(i+2) 


Hence, the following equeation is applied to calculate the fitness function. 

4.3 Best Set formation 

In the proposed HPSO, to improve efficiency of the algorithm, the best solutions which are 
obtained so far are selected and kept in the Best Set. Then, the Best Set is applied to forming 
the Frequency Matrix in next phase of the algorithm. 

To form the Best Set, in the first iteration of algorithm and after the cloning phase of the 
algorithm, the B first best particles among all particles in the swarm are selected and placed 
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in the Best Set. In the other iterations, only the particles that better than the existed particles 
in the Best Set are replaced with them. 

4.4 Frequency Matrix formation 

The Frequency Matrix is a matrix which represents the average times that a specific job goes 
to a specific position according to sequence of particles in the Best Set. To illustrate the 
Frequency Matrix formation procedure, assume that the following particles are in the Best 
Set. 

First particle (sequence): (1,2,3,4,5) 

Second particle (sequence): (1,2,4,3,5) 

Third particle (sequence): (1,2,3,5,4) 

Therefore, the Best Set will be as follows (Figure 6): 


^^MJosition 

Job 

1 

2 

3 

4 

5 

1 

1 

0 

0 

0 

0 

2 

0 

1 

0 

0 

0 

3 

0 

0 

.66 

.33 

0 

4 

0 

0 

.33 

.33 

.33 

5 

0 

0 

0 

.33 

.66 


Fig. 6. The example Frequency Matrix 


4.5 Inversion mutation 

The mutation operator causes a random movement in the search space that result in solution 
diversity. Inversion mutation is adopted in the proposed algorithm. The inversion mutation, 
as illustrated in Figure 7, selects two positions within a chromosome at random and then 
inverts the subsequence between these two positions. 


1 

8 

2 

5 

4 

3 

6 

7 

1 

8 

3 

4 

5 

2 

6 

7 


Fig. 7. General scheme inversion mutation 
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5. Experimental Results 

The performance of the proposed hybrid particle swarm optimization is compared with 
three well-known metaheuristic algorithms: GA, PSO-I, and PSO-II. These algorithms have 
been coded in the Visual Basic 6 and executed on a Pentium 4, 1.7 GHz, and Windows XP 
using 256 MB of RAM. Note that the performance of the proposed algorithm is also 
compared with Lingo 8 for small-sized problems. 

5.1. Benchmark algorithms 

At first, we present a brief discussion about the implementation of benchmark algorithms: 
GA, PSO-I, and PSO-II. 

5.1.1 Genetic algorithm (GA) 

Genetic Algorithm (GA) was developed by Holland in 1975 as a tool for solving complex 
optimization problems of large solution search spaces (Holland, 1992). GAs have been 
applied successfully to a wide variety of optimization problems to find optimal or near- 
optimal solutions (Gen and Cheng, 1997). Thus, for evaluating the performance and 
reliability of the proposed PSO algorithm, we use GA as one of three benchmark algorithms. 
A pseudocode for the applied GA is provided in Figure 8. 

Begin; 

Generate random population of N solutions; 

For each solution: calculate fitness; 

For i= 1 to number of generations (G); 

For y=l to N x CrossovrRate; 

Select two parents randomly; 

Generate an offspring = crossover (Parent 1 and Parent2); 

Calculate the fitness of the offspring; 

If the offspring is better than the worst solution then 
Replace the worst solution by offspring; 

Else generate a new random solution; 

Next; 

Do 

Copy the i th best solution from previous generation to current generation; 

Until population size ( N) is not reached; 

For k=\ to N x Mutation Rate; 

Select one solution randomly; 

Generate a New Solution = mutate {Solution); 

Next; 

Next; 

End. 


Fig. 8. Pseudocode for the Genetic Algorithm 
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5.1.2 PSO-I (Basic algorithm) 

In this section, the structure of PSO-I (basic algorithm) is briefly described. The pseudocode 
of the ap plied PSO-I is provided in Figure 9. 


Initialize the particle population randomly 

Do 

Calculate fitness values of each particle 

Update pBest if the current fitness value is better than pBest 

Determine nBest for each particle: choose the particle with the best 

fitness value of all the neighbors as the nBest 

For each particle 

Calculate particle velocity according to (19) 

Update particle position according to (20) 

While maximum iterations or minimum criteria is not attained 


Fig. 9. Pseudocode for the PSO-I Algorithm (Shi and Eberhart, 1998) 


PSO is initialized with a group of random particles and then search for optima by updating 
each generation. In each iteration, particles are updated by following two best values. The 
first one is the location of the best solution a particle has achieved so far which referred it as 
pBest. Another best value is the location of the best solution in all the population has 
achieved so far. This value is called gBest (Shi and Eberhart, 1998). Equation (19) calculates a 
new velocity for each particle as follows. 

V id = w x V id + Cj x Rand () x (pBest id - x id ) + c 2 x rand () x (nBest nd - x id ) ( 19 ) 
Where Rand () and rand () are two random numbers independently generated. C l and 

C 2 are two learning factors, which control the influence of pBest and nBest on the search 

process. The global exploration and local exploitation abilities of particle swarm are 
balanced by using the inertia weight, W . Particles' velocities are bounded to a maximum 

velocity V max for managing the global exploration ability of PSO (Shi and Eberhart, 1998). 


Equation (20) updates each particle's position ) in the solution hyperspace. 


x id= x id + v id 


(20) 


5.1.3 PSO-II (Constriction algorithm) 

In this section, the structure of PSO-II (constriction algorithm) is expressed in a few words. 
Also the structure of PSO-II is similar to PSO-I (as illustrated in Figure 4), but in PSO-II the 
velocity for each particle is calculated according to equation (21) (Engelbrecht, 2005). 


V id = X [ V id + A ( P BeSt id ~ X id ) + <t>2 ( nBeSt nd ~ X id )] 


(21) 


Where 
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2k 

(22) 

0 = <l> l +<f>2 

(23) 

(f) x = c) x Rand () 

(24) 

</> 2 =c 2 x randQ 

(25) 


Equation (22) is employed by considering the constraints that (j)> 4 and k E [0,1] • By 

employing the constriction approach under above mentioned constraints, convergence of 
the swarm to a stable point is guaranteed. The exploration and exploitation abilities of the 
algorithm are controlled by the parameter of equation (22): k (Engelbrecht, 2005). 


Small-sized problem 

Large-sized problem 

No. Of 
Parts 

Problem 

Number 

Problem 

Condition 

No. Of Parts 

Problem 

Number 

Problem 

Condition 

5 

1 

a i > b t > Cj 

50 

22 

a t > b t > Ci 

2 

AI 

AI 

o* 

23 

a i * Ci > b. 

3 

h , * a t > c t 

24 

b, > a t > Ci 

4 

h , * c t > a t 

25 

bi > Ci > a t 

5 

c i > «, ^ b x 

26 

S) 

IV 

IV 

6 

c i ^ b t > a i 

27 

c t > bi > a t 

7 

Unconditional 

case 

28 

Uncondition 
al case 

10 

8 

a, i > b t > c t 

75 

29 

a i > bi > Ci 

9 

AI 

AI 

o* 

30 

a t > Ci > ^ 

10 

h , * a, ^ c t 

31 

^ > a t > Ci 

11 

h , * c t > a t 

32 

^ > Ci > a. 

12 

c t > a t > b t 

33 

S) 

IV 

IV 
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13 

Ci > b i > a,- 


34 

Ci > bi > a t 


14 

Unconditional 

case 


35 

Uncondition 
al case 

15 

15 

a t > b t > c t 

100 

36 

a i ^ k > Ci 


16 

3 

IV 

IV 


37 

Al 

Al 


17 

b t > a.i > Ci 


38 

^ > a t > c t 


18 

b t > Ci > a,- 


39 

^ > Ci > ai 


19 

IV 

IV 


40 

Ci > d t > ^ 


20 

Ci > bi > a,- 


41 

Ci > ^ > di 


21 

Unconditional 

case 


42 

Uncondition 
al case 


Table 1. Problem inctances 


5.2 Test Problems 

To validate the proposed model and the proposed algorithm, various test problems are 
examined. The experiments are implemented in two folds: first, for small-sized problems, 
the other for large-sized ones. For both of these experiments, the values of £ and 8 are 
equal to 1; the processing time for all parts on the all machine are uniformly generated in 
range [10, 100]. The problem instances are randomly generated as Table 1. 

5.3 Parameters selection 

For tuning the algorithms, extensive experiments were accomplished with different sets of 
parameters. In this section, we only summarize the most significant findings: 

Genetic algorithm 

No of Generation, Population Size, Crossover Rate (Linear order Crossover) and Mutation 
Rate (Inversion Mutation) for the small-sized problems were set to 50, 50, 1.0, and 0.2, 
respectively; and for the large-sized problems were set to 100, 100, 1.0 and 0.2, respectively. 

PSO-I algorithm 

No of Generation, Swarm Size, Learning factors ( C x and C 2 ), and V max for the small-sized 

problems were set to 50, 50, 2, 2, and 3, respectively; and for the large-sized problems were 
set to 100, 100, 2, 2, and 3. The inertia weight for all problem inctances was set to 1.4 that 
linearly decreases to 0.9 in each iteration. 

PSO-II algorithm 
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No of Generation, Swarm Size, Learning factors ( C x and C 2 ), and V max for the small-sized 
problems were set to 50, 50, 2, 2, and 3, respectively; and for the large-sized problems were 
set to 100, 100, 2, 2, and 3. For all problem inctances, k was set to 0.5. 

HPSO algorithm 

No of Generation, Swarm Size, Learning factors ( C x and C 2 ), and V max for the small-sized 

problems were set to 50, 50, 2, 2, and 5, respectively; and for the large-sized problems were 
set to 100, 100, 2, 2, and 5, respectively. Mutation Rate, Best Set size. Clone size, and F for all 
problem inctances were set to 0.1, 7, 5, and 3, respectively. The inertia weight for all 
problem inctances was set to 1.4 that linearly decreases to 0.9 in each iteration. 

5.4 Numerical results 

In this section, the proposed HPSO is applied to the test problems, and its performance is 
compared with above mentioned benchmark algorithms. Each algorithm was executed for 
15 times and the mean results were calculated. The numerical results for various test 
problems are presented in Tables 2 and 3. 
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Table 3. Computational results for large-sized test problems 
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As shown in Tables 2 and 3, the proposed HPSO is superior to the benchmark algorithms in 
the most test problems. 

As illustrated in Tables 2 and 3, the proposed HPSO consumes more computational time 
than the benchmark algorithms. Because of the structure of the proposed HPSO, it can 
search smartly more regions of the search space that results in better solutions. Thus, this 
higher value of computational time is reasonable. 


6. Conclusions 

This chapter developed a new mathematical model for a cyclic multiple-part type three- 

machine robotic cell problem under S 6 robot movement policy that minimizes the cycle 
time. The developed model is based on Petri nets and provides a new method to calculate 
cycle times by considering waiting times. It was proved that calculating cycle time under 

s 6 policy is unary NP-complete. Hence, this chapter proposed a new hybrid particle 
swarm optimization (HPSO) algorithm to tackle the problem. To validate the developed 
model and solution algorithm, various test problems with different sizes were randomly 
generated and the performance of the HPSO was compared with three benchmark 
metaheuristics: Genetic Algorithm, PSO-I (basic Particle Swarm Optimization algorithm), 
and PSO-II (constriction Particle Swarm Optimization algorithm). The numerical results 
showed that the proposed HPSO outperforms the benchmark algorithms in the most 
problems, especially for large-sized problems. 
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1. Introduction 

Swarm optimization, swarm intelligence and swarm robotics are the fields considering a 
group of relatively simple individuals able cooperate to perform complex tasks, in 
decentralized manner. The inspiration is found in the first line within animal societies, such 
as birds, ants and bees. Social insects exhibit successful behavior in performing complex 
tasks on the level of the group, and are able to eliminate noise, errors, failure of swarm 
members. These swarms are robust, able to adapt to constant environmental changes in 
conditions of limited communications among members and lack of global data. In the 
context of swarm optimization, the example of Dorigo's "Ant Colony Optimization " (ACO) 
and Kennedy ad Eberhart "Particle swarm Optimization" (PSO) are most known examples 
of applying swarm-based concepts to development of optimization algorithms able to cope 
with hard optimization problems. These algorithms are justifiably called swarm algorithms, 
because they are run asynchronously and in decentralized manner (Benni, 2004). They also 
mimic the stigmergic (communication by dynamically changing environment) behavior of 
swarm of insects. 

PSO is inspired by flocking behavior of the birds searching for food. Although PSO shares 
many common attributes with the field of Genetic Algorithms (GA), such as stochastic 
nature, population of solution candidates, PSO methods, unlike GA use a kind of 
cooperation between the particles to drive the search process. PSO methods have no 
evolutionary operators like crossover and mutation. Each particle keeps track of its own best 
solution, and the best solution found so far by the swarm. It means that the particles posses 
own and collective memory, and are able to communicate. The difference between the 
global best and personal best is used to direct particles in the search space. 

ACO employs the search process that is inspired by the collective behavior of trail deposit 
and follow-up, which is observed within real ant colonies. A colony of simple agents, the 
ants, communicates indirectly via dynamic modifications of their environment (trails of 
pheromones) and thus proposes solution to a problem, based their collective experience. 
Honey Bees Mating Algorithm (HBMA) can also be observed as a typical swarm based 
approach to optimization. The algorithm is inspired by behavior of eusocial insects, which 
are characterized by three main features: cooperation among adults in brood care and nest 
construction, overlapping of at least two generations, and reproductive division of labour. 
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respectively. In a recent work, Abbas proposed an optimization algorithm based on honey- 
bees mating process (Abbas, 2001; Abbas, 2002). 

The path planning problem of a mobile robot is to find a safe and efficient path for the robot, 
given a start location, a goal location and a set of obstacles distributed in a workspace 
(Latombe, 1991.). The robot can go from the start location to the goal location without 
colliding with any obstacle along the path. In addition to the fundamental problem, we also 
try to find a way to optimize the plan, i.e. to minimize the time required or distance 
travelled (Du et al., 2005; Sadati and Taheri,2002; Ramakrishnan and Zein-Sabatto, 2002). 

The popular methods are the visibility graph algorithm and the artificial potential field 
algorithm. However, the former lacks flexibility and the latter is prone to suffer from 
difficulties with local minima (Alexopoulos and Griffin, 1992; Chen and Liu, 1997). Neural 
network and genetic algorithm have been shown to be very efficient in robot navigation 
(Zarate et al., 2002). General path planning methods based on neural network always 
establish the neural network model for a robot from the start position to the goal position 
and entail much computational time. The input data of the model are the previous distance 
values and position or direction from the sensors. The output data are the next position or 
direction by self-learning process. 

Genetic algorithm is multisearch algorithm based on the principles of natural genetics and 
natural selection (Goldberg, 1989). Genetic algorithm provides a robust search in complex 
spaces and is usually computationally less expensive than other search algorithms. Genetic 
algorithm searches the solution from a population of points and is less likely to be trapped 
in a local optimum. Many results in the literature show the good application of genetic 
algorithm in robot path planning (Khoogar and Parker, 1991; Ram et al., 1994). 

In this chapter, concept of swarm intelligence, as an optimization technique is proposed for 
finding collision free paths in work space containing differently shaped and distributed 
obstacles. Thus, the problem of path planning is considered as an optimization problem, 
whereat collision free paths receive higher fitness values relative to those resulting in 
collision with an obstacle. Performance of HBMA algorithm is compared to the performance 
of a GA developed for the same purpose on two examples, Diophantine equation problem 
and path planning problem. 

Organization of the chapter is as follows: in section 2 we briefly describe colony of Honey 
Bees, as they are in nature. Section 3 describes proposed abstraction and simplification and 
describes core elements of the algorithm. In section 4 and 5 HBMA is compared with GA for 
the first test case, Diophantine equation, and the performances of both algorithms in terms 
of completeness of the solution and speed of the convergence are discussed. In sections 5 
and 6 both algorithms are applied to the second test case, path planning. We conclude with 
section 7 by finally comparing both algorithms and proposing further possibilities of 
improving and testing of the described algorithms. 


2. Structure of a Honey-Bee Colony 

A honey-bee colony typically consists of a single egg laying queen, usually from zero to 
several thousands drones and 10000 to 60000 workers. Drones are the fathers of the colony. 
They are haploid and act to amplify their mother's genome without alteration of their 
genetic composition except through mutation. Workers specialize in brood care and 
sometimes lay eggs. Broods arise from either fertilized or unfertilized eggs, whereby the 
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former represent potential queens or workers, and the latter represent prospective drones. 
The mating process occurs during mating-flights far from the nest. A mating flight starts 
with the dance where the drones follow the queen and mate with her in the air. In a typical 
mating-flight, each queen mates with seven to twenty drones. In each mating, sperm reaches 
the sprematheca and accumulates there to form the genetic pool of the colony. Each time a 
queen lays fertilized eggs, she retrieves at random a mixture of the sperms accumulated in 
the spermatheca to fertilize the egg. 


3. Artificial Model 

The main processes of the algorithm are: mating flight of the queen with the drones, creation 
of new broods by the queen, improvement of the broods by workers, adaptation of workers 
fitness, replacement of the queen with the fitter brood. The mating flight may be considered 
as a set of transitions in a state-space (the environment) where the queen moves between the 
different states in some speed and mates with the drone encountered at each state 
probabilistically, according to (1). 

At the start of the flight, the queen is initialized with some energy content, typically this is a 
random value from range (0,1] and returns to her nest when energy content equals to zero 
or when her spermatheca is full. In developing the algorithm, the functionality of workers is 
restricted to brood care, and therefore, each worker may be represented as a different 
heuristic which acts to improve a set of broods. 

A drone mates with a queen probabilistically according to annealing function: 

A(/) 

prob(Q,D) = e ^ 

Where prob(Q,D) represents the probability of successful mating, i.e. the probability of 
adding drone's D sperm to queen's Q spermatheca. ^(f) is the absolute difference between 
the fitness of the drone and the queen, and S(t) is the speed of the queen at time t. 

According to defined annealing function, the probability of mating is high when either the 
queen is the start of her flight, and therefore, her speed is high, or when the fitness of the 
new potential drone is similar to the queen's fitness. The main steps of the algorithm are 
presented in Fig. 1. 

After each transition in space, the queen's speed S(t) and energy E(t) decay using the 
following equations: 

S(f + l) = a-S(0 
E(t + \) = E{t)-y 

Where a is a factor in range [0.5, 1] and y is calculated according to expression: 

w M 


(2) 

(3) 


And M is the size of spherma theca. 


50 


Swarm Robotics, From Biology to Robotics 


Define input parameters 


Generate queen's 
chromosome at random 


Use worker (heuristic) to increase 
fitness of Queens chromosome 


Generate list of 
drones at random 



Use simulated annealing to fill the 
spermatheca and form the mating pool 




r 


Generate brood by randomly selecting sperm 
from spermatheca and performing haploid 
crossover with queens chromosome 



r 



Improve new solutions applying 
different heuristics (workers) 



Update fitness value of different 
workers 


Replace the queen with the best 
brood 

Yes ^"Best brood better^N 0 
\than the queen 

Keep the previous queen 








Fig. 1. Flowchart of HBMA algorithm 


4. Algorithm Application to Diophantine Equation 

In order to perform initial test of the algorithm, we apply the HBMA to a benchmark 
Diophantine problem. Diophantine equation is an algebraic function (Bull et al., 2006) which 

must be solved over the integers X. eZ. Diophantine problems have a long pedigree in 
number theory. They also constitute some of the hardest problems in modern mathematics. 
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Behavior and results of HBMA and GA applied to the Diophantine nonlinear equation, i.e. 
Markoff equation: 


2 2 2 ^ 

x +y +z =3 xyz 


( 5 ) 


which has important applications in number theory and known solutions. This example is 
chosen because it is known how to generate all the solutions in a cube of given size. In the 
first test case, the problem is reduced to a 2D space by fixing z=433, to have a unique 
solution, and finding integers that satisfy: 

x 2 +y 2 + 433 2 - I299xy = 0 ( 6 ) 

with the search space highly complex in size, as presented with Fig.2. 


F(x,y) 



Fig. 2. Search space for the reduced Diophantine problem 

5. Results for Diophantine Equation 

HBMA and GA were applied to find solutions of the described problem by searching for 
values in the range [0,400]. Both algorithms were successful finding solutions for the 

problem, resulting with monotonous shape of the fitness functions, as presented with Fig. 3. 
and Fig. 4. 
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The fitness function equals the value defined with eq. (6) and is normalized to the range 
[0, l] . In other words, pairs of numbers which yield lower values of fitness function have 

higher chances of survival, ideally approaching zero value for solution and termination 
criteria satisfaction. 



Fig. 3. Fitness value for the HBMA for Diophantine equation 



Fig. 4. Fitness value for the GA for Diophantine equation 

It is important to notice that performance of the HBMA depend on the depth of stochastic 
search. In this example, only two workers i.e. different heuristics were included; namely, 
random walk (RW) and two point crossover (2PCO). That means that in each generation of 
the main loop, a number of local iterations (heuristics) take place for improvement of the 
brood. In our examples, depth of the local searches is set to 100 iterations. 
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HBMA has implicitly included elitist function, because the queen is always represented by 
the best chromosome found so far over all previous generations. 

For the GA, 10% elitism is included, meaning that 10% of best chromosomes are directly 
copied to the new generation, resulting with keeping of best genetic material through the 
whole evolutionary search. Results and behavior of HBMA and GA are presented with 
Fig.3. and Fig.4. Fitness values are normalized, and it is possible to directly compare fitness 
values. In the case of HBMA, the search starts from initial value 0.45, what is likely to be the 
consequence of the first worker applied to initial queen's chromosome. 



No. of runs 

No. of solutions 
found 

Average No 
of generations 


GA 

30 

28 

340 

44.6 

HBMA 

30 

30 

22.9 

4.6 


Table 1. GA vs. HBMA performance comparison 

Both algorithms were tested for 30 runs, with results presented in the Table 1. It could be 
summarized that HBMA outperforms GA in terms of the completeness of the solution. In 
terms of speed of the convergence, one should bear in mind that HBMA has inner loop with 
different heuristics. In each generation, there are number of iterations, defined by the depth 
of stochastic search, taking place. 

Parameters of the GA are: crossover probability: 0.7, mutation probability: 0.01, population 
size: 30, survival selection: generational, initialization: random, termination condition: 
solution found or no fitness improvement over the last 50 generations. 

Parameters of the HBMA are: spherma theca size: M=12, stochastic search depth: 100, 
number of broods: 30, queens energy E and speed S randomly initialized on range [0.5,1], 

0.5 E(t) 

energy reduction step y = , heuristics included: random walk and two point 

M 

crossover. Termination conditions are: solution found or no queens fitness improvement 
over last 50 generations. 

6. Path Planning Results 

HBMA algorithm is implemented to solve the problem of navigation of the mobile robot 
through the space containing arbitrarily distributed obstacles. The environment presentation 
is based on occupancy grid representation. Occupancy grids represent the world as a two- 
dimensional array, with each cell having particular value of 1 (if occupied) or 0 (free cell). In 
our study, obstacles are presented with pairs of nodes connected by mathematically defined 
lines. This is a more compact way of presenting obstacles which will be shown as very 
useful for determining collisions with the KBA. It is possible to create different obstacles as 
lines, or polygons, both convex or concave easily using this compact representation. To be 
able to treat the mobile robot a point in the environment, a minimum safety distance is 
added on the nodes producing a safety shadow around the actual obstacles. 
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Fig. 5. Environment presented in form of occupancy grid. Numbers denote different nodes. 
Bold lines: obstacles; solid lines: feasible path; dashed line: unfeasible path. 


One possible mobile robot environment is presented in Fig. 2. Obstacles are defined as lines 
connecting corresponding nodes e.g. nodes 21 and 23 are occupied and connected with the 
first line, making the intermediate node 22 also occupied. Nodes 55 and 47 are connected 
with the second line etc. Lines can create different shapes, making nodes falling into the 
polygons, " unavailable" for the robot. In case of vertical lines, which cannot be defined as 
mathematical functions since mapping x— >y is not uniform, a threshold value is defined 
such that threshold — ► 0 and added on the x value of second boundary node of the line. In 
such manner, line is slightly rotated around the first node, without real impact on the 
obstacle position and mathematical consistence is preserved. 


6.1 Objective Function 

Impact of the objective or fitness function has a crucial role on the overall performance of the 
evolutionary-based algorithms. The main concept in evolutionary robotics has so far been 
the definition of an effective fitness function (Mermigikis & Petrou, 2006). The authors 
propose some kind of methodology and state that in order to achieve evolution of useful 
behaviours, the corresponding fitness function must have the simplest possible form 
(implicit), it must be possible to be calculated by means of the robot itself (intrinsic) and 
includes elements of the behaviour itself rather than functional details of how this can be 
achieved. Proper form and tuning of the parameters can significantly increase speed of the 
convergence and reduce the possibility of trapping in local optima. In evolutionary-based 
algorithms, objective function has the role of selection of individuals competing to be 
selected for the breeding pool and to transfer their genetic material to the new population 
through the offspring. In the problem being in focus here, the objective function has to 
reward those individuals (paths) that result in minimal number of collisions with obstacles 
and travel minimal distance from the start to goal position at the same time. Fitness function 
is presented by eq. 7: 
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Fitness 


value = - 


* X"/ XV 


7=1 


(7) 


Where wl and w2 are weight constants, w l a w 2 > 0 , X « • is number of collisions of current 
trajectory with obstacles; X<T is total Euclidean distance travelled from origin to 

destination point for current trajectory composed of n components; A is a constant and 
A > 0 . 

Fitness function penalizes trajectories resulting with more collisions and larger total distance 
travelled. To check collisions of the trajectory i and obstacle k, two cases can occur. Case 1: a 
going-through node falls onto the obstacle. This situation is easy to detect and to handle. 
Case 2: a part of the trajectory between two consecutive going-through nodes intersects 
obstacle. This case is handled by solving linear systems of equations for each line segment of 
the trajectory and for each obstacle as a result of following system of presented with Eq. 8. 


X = A 1 B (8) 

Matrices A and B contain coefficients derived from lines that describe obstacles and line 
segments of current path. Matrix X contains solution of linear system and contains point of 
intersection of obstacle and linear segment of the trajectory. If intersection point of any line 
segment S and any obstacle O lies on that particular line segment, then trajectory x intersects 
obstacle O. Otherwise obstacle O doesn't intersect trajectory x. 

Formally: 

© = U o t 

i = 1 

Pi = y> y< (*-*<) 

l x > x < 




y-y<i 

y>i ~y<i 


O t := {(x,_y)|x el;0< (y) < l} 


T 


k 




r[~l©= Ur 00,. 

i=i 

rno, = us 1 , no, 

j = i 

S J ^Pi={xj i ,y jt ) 

Intersection of segment S t and obstacle O t : 

. Vi no 7(r„i' i 

0 otherwise 


(9) 


Values of weight factors are environment dependent and determined experimentally in this 
study, although parameterization of environment with regards on number and distribution 
of the obstacles is considered for future work. This parameterization will include number of 
obstacles, distribution (spread or clustered) and position of obstacles in environment (along 
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the path connecting initial and goal position, or in corner away of main pathways). Through 
parameterization, correlation of form of objective function, neural architecture and 
presented environment could be revealed and thus efficiency of the algorithm further 
increased. 

6.1 Simulation Results 

Different environmental setups were used for the experiments. Performance of both 
algorithms significantly depends on the distribution of the obstacles, namely, whether 
obstacles are cluttered, concentrated, in the vicinity of the goal position etc. The most 
difficult environmental setup is when obstacles are cluttered around the proximity of the 
goal position. 


Environment 



Fig. 6. Environment with obstacles, dashed, unfeasible path, red and feasible path in green 
colour 

One possible environment setup is presented with Fig. 6. Four obstacles are present in the 
environment with given initial and destination position. For the environment presented 
with Fig. 6., comparison of HBMA and GA is conducted. Results are presented in the Table 2. 



No. of runs 

No. of solutions 
found 

Average No 
of generations 


GA 

500 

478 

3120 

430 

HBMA 

500 

493 

430 

58 


Table 2. GA vs. HBMA performance comparison for path planning problem. 

For simplicity, 10 x 10 grid is applied to the environments. Parameters of the GA are: 
Population size = 50, crossover probability: 0.8, adaptive mutation rate: start with 0.1, 
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increment 0.1 if no fitness improvement over 50 consecutive steps. Selection is roulette- 
wheel generational, with the best member of previous generation replacing the worst 
member of current population. Maximum length of chromosomes (degrees of freedom of 
trajectory) =15. 

Both algorithms are able to find solutions for the presented environment with relatively 
high confidence. Again is the completeness (total number of the solutions found by the 
algorithm) slightly on the side of the HBMA. At the same time, number of iterations 
required is lesser for the HBMA, but CPU time is larger, because of the presence of the 
internal loop for the brood improvement. 

Parameters of the HBMA were the same as in the Diophantine equation example. Regarding 
the problem of appropriate parameter selection, it is known to be difficult to tune 
parameters for optimal algorithm behavior, for both algorithms. Parameters were 
experimentally chosen.. 


7. Conclusions 

HBMA algorithm was developed and compared with performance of the GA algorithm for 
two test cases. The firs test case was a benchmark Diophantine equation problem. It is 
shown that HBMA is comparable to the performance of well known GA in terms of CPU 
time, with the time slightly on the side of the GA. In terms of completeness of the solution, 
HBMA was able to find all solutions for the given problem, whereas GA twice did not find 
the solution for given termination criteria. 

Similar behavior was observed for the second test case, namely collision free path planning 
for the mobile robot. However, it is not easy to conclude that HBMA outperforms GA in any 
way, since both algorithms are stochastic and dependant on the proper selection of 
parameters. Although both algorithms and objective were designed to be as simple as 
possible, to enable fair comparison, additional experiments should be performed to achieve 
more reliable behavior and merits for the algorithms. 

HBMA could be further improved by adding additional workers (heuristics) and by 
monitoring success of different heuristics on different problems. GA could be improved by 
tailoring specific evolutionary operators for given problems. 
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1. Introduction 

Extending the particle swarm optimization (PSO) algorithm to be one of systemic modeling 
and controlling tools, several research groups investigate target search with swarm robots 
(simulated or physical) respectively (Doctor et al., 2004; Hereford & Siebold, 2008; Jatmiko 
et al., 2007; Marques et al., 2006; Pugh & Martinoli, 2007; Xue & Zeng, 2008). The common 
idea they hold is to map such swarm robotic search to PSO and deal it with by employing 
the existing bio-inspired approaches to the latter case in a similar way (Xue et al., 2009). Of 
the mapping relations, some aspects including fitness evaluate and path planning have to 
be especially considered because PSO-type algorithm working depends heavily upon them. 
Unlike regarding these respects in PSO, however, the actual characteristics of robot and com- 
plexity of sensing to environment make it impossible to be simplified even ignored. Bear that 
in mind, we might as well explore some representative research work. Pugh et al. compare 
the similarities and differences in properties between real robot and ideal particle, then extend 
PSO directly to model multiple robots for studying at an abstract level the effects of changing 
parameters of the swarm system (Pugh & Martinoli, 2007). Xue et al. simplify characteristics 
of robot by treating each physical robot as a first order inertial element to study mechanism 
of limited sensing and local interactions in swarm robotic search (Xue & Zeng, 2008). Doc- 
tor et al. discuss applying PSO for multiple robot searches, whose focus is on optimizing 
the parameters of their algorithm (Doctor et al., 2004). Jatmiko et al. exert mobile robots for 
plume detection and traversal, with utilizing a modified form of PSO to control the robots 
and consider how the robots respond to search space changes such as turbulence and wind 
changes (Jatmiko et al., 2007). Hereford et al. consider how well the PSO-based robot search 
will scale to large numbers of robots by designing specific communication strategies. Based 
upon this, they have published results of implementing their PSO variants in actual hardware 
robot swarms (Hereford & Siebold, 2008). Marques et al. analytically compare PSO-based 
cooperative search and gradient search as well as biased-random walk search to try to find 
out which performing well in search efficiency. Due to the exchange of information between 
neighbors in the first search mode, PSO-type olfactory guided search possesses merit in search 
properties over its two competitors (Marques et al., 2006). 

It is clear that all of works mentioned above neither involve target search with PSO-type con- 
trol algorithm under conditions of realistic sensing to environment, nor handle the problem 
of obstacle avoidance in the process of target search. On the contrary, each of them assumes 
a potential target in search space to give off a diffuse residue that can be detected by a single 
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sensor, which not corresponding with the actual needs and only having theoretical signifi- 
cance (Hereford & Siebold, 2008). In fact, target signals in the real world can not simply be 
attributed to only one type. Thus, it is need to treat real-time heterogeneous signals fusion 
rather than measure unisource signals as fitness evaluate. Just take search and rescue in disas- 
ters for example. When working miners are confronted with gas outburst accidents in closed 
roadways, they would be likely to lose touch with outside. Unfortunately, search operations 
here tend to be difficult because of the extreme risk. Swarm robots may therefore be pitched 
into carrying out such missions taking the place of human beings. There are multiple kinds 
of heterogeneous signals, including intermittent sound of call for help and periodic radio fre- 
quency (RF) waves as well as continuous gas on disaster spot. We thereupon conduct a case 
study of target search for propose of PSO-type control. Thus, the rest of this paper proceeds as 
follows: In Section 2 the system modeling at individual and swarm levels is done to introduce 
the topics. In Section 3, the properties of target signals are introduced, then a fusion frame- 
work is presented. Then, real time path planning strategy for a typical swarm of wheeled 
mobile robots (WMR) with kinematic constraints in unstructured environment is described 
in Section 4. To examine the validity of fusion approach and path planning, simulations are 
conducted in Section 5. Finally, we conclude in Section 6. 

2. System Modeling 

Consider a swarm of N differentially steered WMRs. The reactive control structure of robot 
used here makes environment sensing linked with actions directly, without requirement for 
explicit expression about search space. Meanwhile, the model of our swarm robotic system 
can be given after mapping PSO to swarm robotic search (Xue & Zeng, 2008). Obviously, the 
modeling to the system consists of two levels according to abstract degrees, i.e., the micro- 
scopic (individual) and the macroscopic (swarm) (Lerman et al., 2005; Martinoli & Easton, 
2002; Martinoli et al., 2004). 



Fig. 1. Control Hiberarchy of Individual Robot 


2.1 Modeling for Individual Robot 

As the reactive architecture with three functional modules including environment sensing, 
behavior planning and actuator driven is chosen (Murphy, 2000), see Fig. 1, both proximity 
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sensor system modeling and kinematic modeling should be considered, while the modeling 
for target signals detection system is dismissed. The reason is that the two former parts are 
related to path planning and the latter configuration depends on the specific types of target 
signals rather than physical size of robot and its kinematics. 

2.1 .1 Proximity Sensor Systemic Model 

To integrate collision avoidance mechanism, we assume that proximity sensors (infrared or 
laser) are equipped on each robot (Jatmiko et al., 2007). Without taking the types and proper- 
ties of proximity sensors into account, we can only extract the commonness according to the 
principle of range measurement for modeling. As for the specific configuration of proximity 
sensors in this work, sixteen proximity sensors are assumed to be equipped on each individual 
robot, surrounding their body in a discrete circular uniform distributional fashion, see Fig. 2 
and Tab. 1 for details, where the black ovals stand for proximity sensors and the arrow points 
the heading. 



Fig. 2. Proximity Sensor System of Individual Robot 


2.1.2 Kinematic Model 

The modeling of the kinematics of robots in a two-dimensional plane can be done using either 
cartesian or polar coordinates. The model in cartesian coordinates is the most widely used 
and discussion here will be limited to modeling in cartesian coordinates (Maalouf et al., 2006). 
Typically, the posture of robot at any instant is defined by the position and heading relative to 
the global frame. The kinematic model is given as follows (Campion et al., 1996): 

{ ±i = Vi cos 6j 

Vi = V{ sin $i (1) 

6i = 

where p; = (xj, y ; ) T be position vector or cartesian coordinates of robot R, under global frame, 
Qj its orientation or steering angle, Vj translational or driving or linear velocity and a ; z - angular 
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Table 1. Distribution Degrees of Proximity Sensors 


62 


Swarm Robotics, From Biology to Robotics 


or steering velocity, as is shown in Fig. 3. In the absence of obstacles, the basic motion tasks 
assigned to a WMR may be reduced to moving between two robot postures and following a 
given trajectory (Oriolo et al., 2002). Whichever can finally be attributed to the design of con- 
trol laws, i.e., control command series of inputs (z;,a;) T . Although the actual commands may 
come in another forms, e.g., the angular velocities cor and coi of the right and left wheels, 
respectively, rather than v and to, we can still make use of analytical module built in robot 
controller to get the required commands by a one-to-one mapping between these velocities 
(Oriolo et al., 2002; Siegwart & Nourbakhsh, 2004). For all control schemes, in fact, an addi- 
tional filtering of original velocity commands is included to account for robot and actuator 
dynamics. 



X 


Fig. 3. Kinematic Model of Individual Robot 

Due to the abilities of motion mechanism and actuators, there exist limitations on real veloci- 
ties. Consequently, the actual input commands can be obtained with the following rules: 


Vmax/ 

if Vj(t) > v max 

0 , 

t-K 

A 

o 

Vi(t), 

otherwise 

C^max/ 

if (Vi(t) > CO max 

~ CO max/ 

if U>i(t) < —COmax 


otherwise 


Note that above rules come from non-holonomic constraints because robot can move along its 
bearing only, that is, the direction of V{ is always in accordance with the heading of robot, see 
Fig. 3 and Fig. 4. Then, can be used to decide the orientation of robot Rj. As for the robot at 
position pi with V((t) = (pn, u# )tr we can calculate the orientation: 

6 At) = arctan Vl2 ^\ 

(0 


( 3 ) 
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Fig. 4. Motion Control of A WMR 


Similarly, the orientation 6{(t + At) = f oc at position p 2 with + At) = (^z'I/^Wa* can be 
decided too. Therefore, the required expected turning angle from p\ to p 2 can be computed 
and used to further decide cop 


A fa = arc tan 


0(2 (t + At) 
v a (t + At) 


— arctan 


vn(t) 


( 4 ) 


The posture vectors of robot R, are required as control inputs to individual con- 

troller at each time step, depending on the posture estimate with incremental encoder data 
(odometry). Assume the angular wheel displacements having been measured during the sam- 
pling time At by the encoders. We can further obtain the linear and angular displacements As 
and AO. Then, the estimate of posture at time t + At can be computationally decided (Oriolo 
et al., 2002; Siegwart & Nourbakhsh, 2004): 


/ x \ / x \ / Ascos(0+^) 

y = y + Assin(e + f) 

V ^ / t+At V (> ) t V Ae 

2.2 Modeling for Swarm Robots 

Our swarm robotic system is composed of the above-mentioned robots. The meanings of 
used symbols are as follows: pi = (xii/X^) and = (vn , vq) are position and linear velocity 
of robot Rj at time t respectively; p* = (x* l7 x* 2 ) and p^ = the best historical 

positions of robot i itself and its communication-based neighborhood (Pugh et al., 2006). Based 
on this, we can define the best position within its neighborhood (Xue & Zeng, 2008; Xue et al., 
2009): 



P*,-)(0 = P*V)/ ar gk max {f(PifcW)/ k € R i' s neighborhood} (6) 

where I() is the fusion of measurement readings of target signals. Further, we are able to 
model swarm robotic system with the extended PSO method: 
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+ 1) = £iVij(t) + c\r\{x*- - Xij) + c 2 r 2 (x* i); . - x;y) 
Vij(t + At) = Vjj(t) + Ki(vjj(t + 1) — Vij{t)) 

Xij{t ~\~ At) — Xij(t) Vij(t At) At 


( 7 ) 


where i? z -y (£) and x z y (£) are /-dimensional velocity and position of robot Rj at time t respectively, 
Vjj(t + 1) the expected computational velocity, K z designed parameter of local controller gain 
which can be chosen by designer. As robot may have to take several time steps AA t, in most 
cases, to reach an expected position from the consecutive previous expected one, adding this 
factor to obtain a "smoother" displacement. Besides, £ z be algorithmic inertia coefficient that 
can be set constantly or dynamically, C\ and c 2 cognitive and social acceleration constant re- 
spectively, r\ and T 2 stochastic variables subject to the distribution of 11(0, 1). We can calculate 
the linear velocity: 



(8) 


With kinematics model given in Eq. (1), the input of linear velocity (Peng & Akella, 2005) and 
angular velocity (Siegwart & Nourbakhsh, 2004) can be decided: 


Vj (t + At) 
co 1 (t - \- At) 



( 9 ) 


'max 


where /3 Z is the computational expected turning angle from the current position to the next 
expected one, also see Fig. 4 for details. 

3. Signals Fusion 

The key to PSO-type search algorithm is to take detection and fusion of target signals as fitness 
evaluate so as to decide the best-found position, since each individual robot is guided by 
the best experience of itself own and its neighborhood. Obviously, the temporal and spatial 
features of different types of signals should be explored in advance. 

3.1 Signal Properties 

With mathematical models of signals propagation, we can generate a set of theoretically com- 
puted signal strength data akin to the empirical data set to design fusion algorithm rather than 
collect the data on spot. 

3.1.1 Sound 

The identical model may be used for both propagation and measurement as to the same space. 
Compared with the size of environment, the mouth of victim can be considered as a point 
sound source. Let N robots equipped with acoustic sensors construct a mobile sensor field, 
where an immovable target emits omnidirectional acoustic signals. The signal energy mea- 
sured on the i ih sensor over time interval t, denoted by (Li & Hu, 2003): 



(10) 
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where f z is time delay for sound propagates from target to the z th robot, s(t) is a scalar denoting 
energy emitted during sampling time t; r(t) coordinates of target during t; r z - coordinates of 
the z th stationary sensor; gj gain factor of the z th acoustic sensor; a(« 2) energy decay factor, 
and £i(t) cumulative effects of modeling error of gj, r z , oc and the additive observation noise 
of y z (t), see Fig. 5. 



Fig. 5. Acoustic Energy Loss Fitting 


3.1.2 RF Waves 

Typically, underground mine personnel tracking systems rely more and more upon radio fre- 
quency identification (RFID) technologies today. Such a system has basic components includ- 
ing readers and tags. The latter is categorized as either passive or active (Ni et al., 2004). As 
for coal mine application, a tag is often mounted on a miner's helmet with his lamp. For a ra- 
dio channel, the transmitted signal reaches receiver via multiple paths (Bahl & Padmanabhan, 
2000 ): 


p(d) = p(d 0 )-io*ig4--n (ii) 

where oc indicates loss rate, P(do) is signal power at reference distance do and d transmitter- 
receiver distance. The value of P(do) can be derived empirically or obtained from the wireless 
network hardware specifications. In general, rj value is derived empirically, see Fig. 6. 

3.1.3 Gas 

The gas in coal mines will diffuse quickly in closed roadways after gas outburst. The pervasion 
process can be described as affecting by some odor point sources. For convenience, they can 
be viewed as only one by linear combination. Let the projection of leak point on the ground 
be origin, average direction of downwind x-axis, a right hand three-dimensional coordinate 
system can be set. Then, we can calculate gas concentration in any point on the ground (z = 0) 
following the law (Marques et al., 2006): 

Q exD r (y(0-yo(*>0 ) 2 1 

2na y (x,t)o- z (x) -2 o$(x,t) 1 


C(x,y,t) = 


(12) 
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Fig. 6. Log-Normal Signal Loss Distribution 


where Q represents release rate, plume center yo, width w and height h as a function of time 
t and downwind distance x. dy{x,t) = , Pz(x) — Fig- 7 shows an example of 

a time averaged Gaussian plume (Marques et al., 2006). Further, the heterogeneous signals 
distribution in search environment can be shown in Fig. 8. 



Fig. 7. Gas Concentration Contours on the Ground after Emitted Sufficient Long Time 


3.1.4 Signals Propagation Space 

The space is divided into six sub-areas (numbered Area 1-6) according to the distribution of 
signals. The lines in Fig. 8 represent the minimum detectable signal contours corresponding 
to thresholds 0.0016 kg/ m 3 , —90 dBm and to maximum detectable ranges 200 m, 45 m respec- 
tively Li & Hu (2003); Ni et al. (2004). It is need to point out that the sound threshold is not 
given definitely because it is closely related to the sensor sensitivity. Hence, given a specific 
power (milliwatt magnitude) of call-for-help in a loud voice, our attention lines in finding 
how far the emitted sound signals can reach. In simulation, we will make an experiential but 
reasonable assumption on this value. 
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Fig. 8. Signals Distribution in Search Space 


3.2 Fusion 

Different types of signals are presented in different forms. For instance, gas diffusion distance 
may be up to several hundreds of meters (Marques et al., 2006); the detectable range of RF 
waves with frequency f = 7.5 s emitted by active tags approach to 150 ft (« 45 m) (Ni et al., 
2004) and the localization accuracy with RF RSSI-based method can be 2 m; the detectable 
range of sound usually reaches not more than 30 m (Li & Hu, 2003), and the estimation error 
with sound RSSI method may be up to 50%. Therefore, detectable range, statistical properties, 
localization types and accuracy have to be considered simultaneously in signals fusion. 

3.2.1 Sensing Event 

Introduce a 2-value logic into describing perceptual process so that we can define perception 
event in advance. Let Aj (i = GAS, RF, CALL) be perceptual event in event space O = {0, 1}. 
Then Aj = 1 represents detect-success (beyond threshold) and Aj = 0 detect-failure. At the 
same time, we normalize those measurement readings beyond threshold to Nm j <E (0,1). 
Clearly, events Aj and Aj (i j ^ j) are mutually independent and the probability of each event 
can be calculated with its statistical properties. Further, let (A zGAS (t), A zRF (t), A zCALL (t)) be 
the joint sensing event of robot Rj at time t, then there may be 2 3 = 8 joint events according to 
the time characteristic of signals distribution. Again, consider the spatial distribution of sig- 
nals. Those robots in signal blind area (Area 1) attempt to capture signal clues independently 
in a spiral move manner to further search for target locally (Hayes, 2002; Marques et al., 2006), 
without directed by swarm intelligence principle locally, while the robots in Area 2-6 do so 
(Marques et al., 2006). Thus, it is easy to know that there are six joint events everywhere ex- 
cept source. The possible joint events occurred in each sub-area are listed in Tab. 2. These joint 
events can be encoded with 3-bit binary numbers. Since the source is characteristic of such 
encode, we can also express it with 3x1 "characteristic" vector C. Finally, take the above 
Nm,s of measurement readings to replace the corresponding elements "1" in vector V. 
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3.2.2 Virtual Communication 

View signals measurement as continuous communication between robots and target. In this 
case, each robot possesses its own channel, with source (target) and destination (robot). The 
individual robot discretizes three types of measurement readings to 1-bit binary digit respec- 
tively with corresponding threshold. As for the same signal emitted from source, robots in 
different sub-areas may obtain different results due to the effect of distance. Consequently, 
while there are 2 3 = 8 encodes of full permutation in source, the " received" encodes by robots 
in different sub-areas may vary. We suppose the duration of emission by target can guarantee 
the detection success in sampling period (400 ms here). But the transmission time from source 
to destination is small enough so as to be ignored. Thus, the detection process can transfer 
continuous information to time- and amplitude-discrete random signal series. Since we sup- 
pose At be sufficient small interval, the above perception event occurs once at most in every 
sampling process. 

3.2.3 Information Entropy 

We can calculate the information entropy from the received information encodes through vir- 
tual communication process. 

• Gases. Robots start to locally search for target as soon as gas can be sensed in 
global search stage (Marques et al., 2006). As to the continuous gas diffusion, event 
Aqas — 1 vvill occur at any time t in Area 2-4 (Hayes, 2002; Marques et al., 2006). Then 
P{ Aqas (f) } = 1/ be., this information is decisive, say, H(Xgas) — 0. 

• RFID Waves. The perception process of RF signals (X RF (t),t > 0} can be viewed as 
Poisson process with intensity Arp. As the process has stationary independent incre- 
ment, those equal intervals may lead to equal probability of A RF = 1. Suppose the 
sampling period At is sufficient small time interval and satisfies the sampling theory. 
The event A RF = 1 occurs once at most in every sampling. Since events A RF = 1 
and A rf = 0 are contrary ones, then we can draw a conclusion P{A RF (t) = 0} = 
1 — P{A RF (t) = 1}. And the relationship can further be captured by computing proba- 
bilities: 


P{A N (t) = !} =e~^ At \ RF At 

p{A RF (f) = 0} = e - A RF Af (13) 

e A RF Af = ;v Rp A f _|_ i 


"Emitted" 

Sub-Area 

Possible Event 

"Received" 

111 

1 

(0,0,0) 

000 

111 

2 

(1,0,0) 

100 

111 

3 

(1,0,0), (1, 1,0) 

100, 101 

111 

4 

(1,0,0), (1,1,0), (1,0,1), 

(1,1,1) 100,110,101,111 

111 

5 

(0,0,0), (0, 1,0) 

000, 010 

111 

6 

(0,0,0), (0,0,1), (0,1,0), 

(0,1,1) 000,001,010,011 


Table 2. Joint Event Encodes (AgaS/ ^RF/ ^call) in Search Space 
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Signal 

Entropy 

Detectable Range(m) 

Localization Type 

Accuracy (m) 

GAS 

0 

200 

indirect 

200 

RF 

0.0156 

45 

direct 

2 

CALL 

0.2055 

30 

direct 

15 


Table 3. Characteristics of Signals Propagated in Search Space 


Accordingly, the entropy of RF signals in information source at t can be calculated with 
the above conclusion: 

H(X rf ) = A RF At + (e ARpAt — 1) log(A RF Af) (14) 


• Sound of Call-for-Help. The detected sound of call for help (X c (f), t > 0} be Poisson 
process with intensity Ac- Similarly, we can determine the probability detect-success 
and detect-failure of call for help to further obtain the entropy of such sound at any 
time t: 


H(X C ) = A C A t + (e- AcAt - 1) log(A c A£) (15) 


3.2.4 Weight 

A criterion to signals fusion can be used by robot. Among all factors, apart from information 
entropy, the localization type and accuracy with RSSI-based method should also be consid- 
ered. For instance, gas source is not same as the target location, i.e., we localize target in- 
directly by localizing gas source based on the fact that one should move quickly upwind in 
risk avoiding poison gas leakage. Consequently, such estimate may get the worst accuracy 
(200 m assumed). While the RSSI-based estimate with RF or sound intensity can localize tar- 
get directly. But the two types of signals differ in accuracy of location estimate, see Tab. 3 for 
details (Li & Hu, 2003; Marques et al., 2006; Ni et al., 2004). As shown in Eq. (16), the entropy, 
localization type and accuracy are all required to be integrated by weighted sums. 

w, = + bK + — i = GAS, RF, CALL (16) 

hH{Xi) T ;Ef 

i i 1 

where logic variable k represents localization type, "indirect" is assigned 0 and "direct" 1. t is 
localization accuracy and a, b, c G (0, 1] are all positive coefficients that need to be determined 
empirically. Then, we can take three weight values as elements to construct a 1 x 3 vector W. 

3.2.5 Signals Fusion 

An fusion mechanism for making decision on the best positions is discussed here, being suit- 
able for deciding on cognitive of individual and social of swarm. The mechanism can be 
expressed with weighted sums operation using vectors V and W, i.e., obtaining the fusion by 
calculating the inner product of two vectors fusion = V(i X 3 ) * W( 3 X i)- 
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3.2.6 Description of Fusion Algorithm 

A full-distributed fusion algorithm run on each individual robot can be presented. First, we 
assume each robot has an unique ID, carrying a set of sensors and a on-board fusion module 
so as to measure signals and fuse them independently. Besides, all sensors are assumed to 
react to signals in sufficient short time. Finally, we design a character structure denoting as 
// ID ,/ +"Position" +" fusion" , which can be viewed as the communication protocol. As for the 
local communication hardware, it can be achieved by wireless transmission systems, like RF 
or infrared. 


Algorithm 1 Real-Time Heterogeneous Signals Fusion 
1: Input: sensor readings 

2: Output: the best-found position and fusion of signals 
3: confirm ID and current position iPos; 

4: initialize 

5: set counter t <— 0; 

6: set (A zGAS , A iR p, A z - call )* =0 = 000; 

7: set fusion = 0; 

8: construct "ID"+"Position"+"/wszoft";{communication protocol} 

9: set best position of itself ibPos iPos ; 

10: set best position of neighborhood sbPos iPos ; 

11: repeat 

12: make measurement; 

13: discretize to 0 or 1 by comparing threshold value; 

14: format data with characteristic structure; 

15: if (A zGAS , A zRF , A zCALL )f = 000 then 

16: keep silence;{do nothing} 

17: else 

18: elect ibPos and update; 

19: broadcast data within its neighborhood; 

20: end if 

21: listen for others; 

22: if receive data containing (Ay GAF , Ay RF , Ay GAFF ) 7^ 000 then 

23: elect sbPos and update; 

24: end if 

25: f<-f + 1; 

26: until termination condition is met 


4. Path Planning 

In PSO-type swarm robotic search algorithm, each individual robot makes decision on its ex- 
pected destination at every time step as its current target to move towards in a full distributed 
fashion by combining its own inertia and cognitive experience as well as experience of swarm. 
The experience of robot itself and its neighbors depends on fitness evaluate, i.e., target signals 
measurement and fusion, which is discussed in the above section. In other words, the trajec- 
tory of each robot "searching" for target is formed by linking a series of expected positions 
orderly. 
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Similar to single autonomous robot, path planning of individual robots in swarm robotic sys- 
tem also involves how to move towards their own goals with static/ dynamic obstacle avoid- 
ance (Warren, 1990). Performed in an iterative manner, artificial potential field (APF) method 
is usually employed for such task in controlling autonomous robot (Khatib, 1986). Theoret- 
ically, robot moves in the direction of the resultant of the attraction force pulling the robot 
towards the goal, and the repulsive force pushing the robot away from the obstacles. As 
expected, the robot stops moving after reaching the goal position. Unfortunately, it always 
suffers from local minima where if trapped (Zou & Zhu, 2003). Looking for a local-minimum- 
free solution has become a central concern in this approach. Aiming at the problem solving, 
some modified APF methods are proposed to overcome local minimum. The key ideas may 
be fallen into two categories: one establishes new potential functions with a few or even no 
local minima (Ge & Cui, 2000; Warren, 1990); the other uses some techniques to escape from 
local minima, including random walk (Janabi-Sharifi & Vinke, 1993), wall following (Boren- 
stein & Koren, 1989), and other heuristic methods (Singh et al., 1996). Except for these efforts, 
Jugh et al. apply the PSO algorithm to path optimization of multiple robots (Pugh & Marti- 
noli, 2006). However, above methods are incompatible with the case of swarm robotic search. 
New difficulties arise when we apply APF method to swarm robots path planning. One major 
challenge is to bridge the high-level task planning and the low-level path planning and inte- 
grate them into one framework (Ren, 2005). Thus we combine APF method and PSO to plan 
path towards target with collision avoidance because of low computational cost and better 
real-time performance. 

4.1 Traditional APF 

The nature of APF method lines in defining the motion space for robot as a virtual potential 
filed U(x) including virtual gravitational field U G (x) and repulsion field U R (x), in which 
robot is attracted by target and repelled by obstacles. Then, the resultant force field can be 
defined with (Khatib, 1986): 


U(x) = U G (x) + U R (x) (17) 

Meanwhile, we can further define attractive force Fq (x) and repulsive force F R (x) as the neg- 
ative gradient of the virtual gravitational field and repulsion field respectively. Therefore, the 
virtual force F(x) acted by the virtual potential field can be derived using space dynamics 
equation and Lagrange equation: 


r F(x)=F G (x)+F R (x) 

< F g (x) = — V (U G (x)) (18) 

{ F r (x) = -V(U R (x)) 

Clearly, the direction of robot motion depends upon the direction of F(x) (Khatib, 1986). 
Although the traditional APF method has the virtue of being the easiest to implement, it 
has some limitations above yet. At first, the traditional APF method is applied to the case 
of global environment information being known rather than the case of environment being 
partly known or even unknown because the virtual potential field is computationally ob- 
tained and what robot sensing environment with its own equipped sensors is not supported. 
Second, the inherent disadvantage of traditional APF method comes through in being easily 
trapped to local minima and in target being not able to reach. It is important that the tra- 
ditional method has to be modified in accordance with robot sensing environment with its 
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sensors. The robot navigates in search space without obstacle collision depends completely 
upon equipped sensors through collecting measurement readings to judge states including 
obstacles distribution and possible target position. 

4.2 Sensor-Based APF 

Generally, when searching for target in unknown environment, the environment map is partly 
known or even unknown. In this case, the robot behaviors for obstacle avoidance have to rely 
on continuous local path planning by means of locally sensing surroundings with equipped 
sensors. As robot moves within search space, the obstacles surrounding robot are inevitably in 
different conditions. Learning from the traditional APF method to improve real-time property, 
we can integrate it with the multi-sensor structure of robot to construct virtual potential force 
with change of sensor readings. Hence, it is need to make some modifications to Eq. (18) based 
on above structural sensor model, see Fig. 2. 



ASij = Sr - Sjj 


where F-(x) be the resultant force imposed on robot Rj in constructed virtual potential field, 
F' g (x) the force attracted by the expected target position, and F' Q (x) the force repelled by 
surrounding obstacles. Furthermore, Sr be the maximum detection range of all sensors and 
Sjj the current distance reading of sensor j, A Sjj represents the increment of the ; th sensor 
reading. Note that A Sjj be a vector because of the directionality of sensors. 

4.3 Control System Architecture 

To decide input commands ( Vj,C 0 j) T of individual robots every time step, the control archi- 
tecture including swarm and individual levels should be deterministic. From swarm aspect, 
the architecture is distributed and the PSO-type algorithm runs on each robot. In individual's 
eyes, robot has a two-level virtual control architecture, which may refer to (Oriolo et al., 2002) 
for details. Our designed algorithm is at high-level layer, running with a sampling time of 
A t = 100 ms. While the low-level layer is charge of analyzing and executing the velocity com- 
mands from upper level. The outputs of algorithm are the command series ( Vj,C 0 j) T in every 
time step. As is shown in Eq. (9), Vj(t + At) and C 0 j(t + At) are the required control inputs of 
linear and angular velocity at the next time step respectively. While Vj(t) and (Oj{t) are the 
obtained current variables by sampling. 

4.4 Description of Control Algorithm 

It is shown that the PSO-type algorithm is capable of controlling individual robots to move 
about in space for target search with obstacle avoidance according to the modified sensor- 
based APF method. Under the conditions of limited sense and local interaction in unknown 
environment, a valid navigation algorithm can be designed for target search with collision 
avoidance. Such idea can be implemented in accordance with the three phases below: 

• Compute the Expected Positions. In terms of the model of swarm robotic system, i.e., 
Eq. (7), the respected velocities and positions of each robot at time step t can be compu- 
tational decided by means of interactions within its own neighborhood. 
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Target 

A 



Fig. 9. Schematic of Virtual Force Acted on Robot with Proximity Sensor Readings 


• Decide Virtual Force. With the modified sensor-based APF model, we can construct 
a potential field and get the virtual force in this field. The specific way is to take the 
expected position of robot at time step as the current temporary target which will attract 
the robot, while the robot will be repelled by the detected static or dynamic obstacles. 

• Compute the Real Positions. As the velocity of robot at time t + 1 is gotten, the position 
of robot at time t + 1 can be computationally obtained according to the kinematics of 
robot. 

A full distributed PSO-type algorithm for target search is developed, which can be imple- 
mented on each robot in parallel. Without loss of generality, we can describe the algorithm 
run on robot R z - as Algorithm 2. 

5. Simulation and Discussions 

To elaborate how to fuse the specific heterogeneous signals and how to decide the best posi- 
tions, the simulations are designed and conducted for the purpose. First, virtual signal gen- 
erators are arranged where same as target situates, emitting signals following their own time 
characteristic. Then, a series of detection points are set in signal Area 1-6. Our task is to inves- 
tigate what happened in each information sink (robot) when different combination of signals 
is emitted from source by virtually measuring and fusing. We observe for sufficient long time 
until all eight encodes transmitted from source. Then we try to find the relationship between 
distance and fusion result. 

5.1 Signals Generating 

Consider the properties of a given Poisson process with intensity A. The successive coming 
time of events obey exponential distribution with mean j . We can empirically set the value in 
some interval, for example, the upper bound and lower bound can set to 0.01 and 0.001 respec- 
tively, i.e., Ac £ (0.001,0.01), while the intensity of RF signals can be Arf = 0.1333 according 
to its primitive definition, which reflect the temporal characteristics of target signals. 
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Algorithm 2 Path Planning for Swarm Robots in A Full-Distributed Way 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 


initialize 

set counter k *— 0; 
initialize constants; 
initialize v l k , x l k ; 
initialize position of target; 
initialize robot's own cognition 
make measurement I l k , 
n <_ n • 

2 max 



initialize shared information 

tS ji. 

%ax ^ 



confirm index of best individual; 

repeat 


fc<-fc + l; 

communicate among neighborhood 
confirm neighborhood; 
for ; = 1 to number _of_neighb or s do 
compute I ] k ; 
fmax 4 max(l|, Jj[); 

pi x r' ar g m max { j ( x r)' m 6 (*'/)}' 

endfor 

compute expected velocity and position 

4+i w fc4 + Wpi - 4) + c 2^(pf - 4); 
4+a^4+ k »(4+i-4)' 

4+Afc^4+4+M M ; 

c 3 ^;{0 < c 3 < 1 } 

compute velocity with kinematics 

4+ai m in(p maX / 4 +ai); 

compute a;| +Afc ; 

if shared information updated by neighbor then 
compute next expected position; 


endif 

until succeed in search 
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5.2 Deployment of Measuring Points 

We set a series of measuring points, assigning one with each sub-area. Different points are 
different far away from the source. Note that a pair of points in different areas having the 
same distance value are arranged to study the relation between fusions at the same time. 

5.3 Main Parameter Settings 

We simulate signals fusion using parameter configuration a = 1, b = 0.001, c = 1, Ac = 
0.2055, Arf = 0.0156. For convenience, target is fixed to (0, 0) all time and the coordinates of 
six measuring points are (100,40), (150,0), (40,0), (20,0), (0,35), (0,20) orderly. Mean- 
while, we focus on if the coverage of all joint events occur in sufficient long time rather than 
the moments. 

5.4 Map Processing 

In study on path planning of autonomous robotics, how to represent the working space, i.e., 
how to model the space is one of the important problems. Based on the difference of sensing 
to environment, modeling approaches to known or unknown map fall into two ones. Here we 
model working space for swarm robots with digit image processing technology. The obstacle 
information relative to each point in search space is expressed with a two-dimensional arrays. 
Of representative symbols, 0 represents passable point and 1 passless. The Fig. 10 is the 
example of mapping processing. 
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(a) Original Map (b) Digitizing 

Fig. 10. Working Space for Swarm Robotic Search 


5.5 Obstacle Avoidance Planning 

Based on the fusion method, we run the swarm robotic search algorithm having a specific 
function of path planning. The unequal sized swarms (N = 3, 5, 8, 10) are used, repeated the 


76 


Swarm Robotics, From Biology to Robotics 


algorithm running for ten times respectively. Then, the statistics about total distance and time 
elapsed in different cases are collected to support our presented method. 

5.6 Results and Discussions 

Conducting the above simulations repeatedly, we can get the following results. And then we 
may hold discussions and draw some conclusions. 

• The fused values in simulation are shown in Fig. 11, from which robots can "find" out 
the best positions by simple election operation. It's perceptible that the bigger the fusion 
value, the nearer the measuring point from target. At the same time, it is observed 
that as for No. 4 and No. 6 points, the fusion results are the same in cases of Source = 
001, 010, Oil, and different in cases of Source = 101, 110, 111 although they are equal to 
distance of target. We may explain it in this manner: robots searching for target depend 
on measurements because they do not know the position of target. While the two points 
are located in different sub-areas, the situation of signals cover is different. 


Source="000" "001" "010" "011" 



" 100 " " 101 " " 110 " " 111 " 



Fig. 11. Fusion results at the assigned six measuring points under different encodes of in- 
formation source. Note that the title Source="000" of the left corner sub-figure represents no 
GAS, no RF, and no CALL signals are emitted when sampling. One can understand the other 
cases in a similar manner. Besides, the fusion is a scalar value without any physical meaning. 


Fig. 12 shows the scenario of two robots decide their respective motion behaviors with 
modified APF method to plan paths for obstacle avoidance. 
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(a) 2-Rob (b) 4-Rob 

Fig. 12. Obstacle Avoidance between Unequal Sized Robots with Sensor-Based APF Method 


Swarm Size 

Average Time 

Average Total Distance 

3 

278 

1930 

5 

232 

2410 

8 

197 

3136 

10 

184 

4380 


Table 4. Statistics from Search for Target Simulations 


• Fig. 13 shows the scenario of one single robot planning its path using multiple sensor- 
based APF method without obstacle collision to search for target successfully under 
different conditions of obstacle types. 

• Consider the total displacements and time (iterative generations) when the search suc- 
ceeds. The statistical results shown in Tab. 4 and the relations between average dis- 
tance/ generations and swarm size are charted as Fig. 14. 

6. Conclusions 

As for PSO-type control of swarm robots, the experience both of individual robots and of 
population is required. In order to decide the best positions, we take the characteristic infor- 
mation of target, such as intensity or concentration of different signals emitted by target, as 
the "fitness". Therefore, the problem of multi-source signals fusion is proposed. To this end, 
we model the process of signals measurement with robot sensors as virtual communication. 
Then, the detected target signals can be viewed as transmitted encodes with respect to infor- 
mation source. We thereupon present some concepts of binary logic and perceptual event to 
describe the "communication" between target and robots. Besides, we also put forward in- 
formation entropy-based fusion criteria and priority to fuse signals and election mechanism 
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(a) Circular 


(b) Rectangle 



(c) Allotype 


(d) Complicated 


Fig. 13. Single Robot Move to the Potential Target with Path Planning 



Fig. 14. Relations between Average Distance/ Generations and Swarm Size 
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to decide the best positions on the basis of space-time distribution properties of target and 
robots. Simulation conducted in closed signal propagation environment indicates the approx- 
imate relation between fusion and distance, i.e., the nearer the robot is far away from target, 
the higher the fusion of signals. Also, a modified artificial potential field method is proposed 
based on the multiple sensor structure for the space resource conflict resolution. The simu- 
lation results show the validity of our sensor-based APF method in the process of search for 
potential target. 
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1. Introduction 

Digital filtering plays an important role in sensors' signal processing of robots. Not like 
analog system, it is not limited by parameters of electronic components, so it can process 
signals of rather low frequency, which is one of its advantages. According to different 
structure, digital filters can be divided into finite impulse response (FIR) digital filters and 
infinite impulse response (HR) digital filters. The output of FIR digital filters only relates 
with the previous and the present input. Whereas the output of HR digital filters relates not 
only with the input but also the previous output. It is to say that HR digital filters have their 
feedback. Seen from signal-processing, HR digital filters have great advantages over FIR 
digital filters, but they also have their disadvantages at design. The coefficient of HR digital 
filters is highly nonlinear, whereas the coefficient of FIR digital filters is linear. 


2. Signal processing system of the robot joint force/position sensor 
2.1 Configure of the signal processing system 

There are two kinds of design methods for HR digital filters: 1) Frequency translation 
method, this method has two design routes: one route is first get analog lowpass filter, 
analog highpass filter, analog bandpass filter and analog band elimination filter by doing 
frequency band transform to the analog normalized prototype, and then get digital lowpass 
filter, digital highpass filter, digital bandpass filter and digital band elimination filter by 
digitization; the other route is first get digital lowpass filter by digitizing the analog 
normalized prototype, and then get digital highpass filter, digital bandpass filter and digital 
band elimination filter by frequency band transform in digital domain. 2) Optimization 
algorithm, it is to design digital filters under certain optimization criterions to get the best 
performance. Now, there are minimum P-error method, least mean square error (LMSE) 
method, linear programming method and model-fitting frequency response method etc. 

In recent years, some scholars have already applied such intelligent algorithms as genetic 
algorithm, artificial immune algorithm and particle swarm optimization (PSO) algorithm etc 
into the design of HR digital filters and achieved better result. Commonly speaking, filters' 
capacity is often shown by the permissible error of amplitude characteristic of its frequency 
response. When designing a filter, we should consider such main technical index as 


82 


Swarm Robotics, From Biology to Robotics 


passband cutoff frequency CO c , stopband cutoff frequency CO c , passband toleranceal, 
stopband tolerance a , 2 and passband maximum ripple , stopband minimum attenuation 

^ 2 , etc. At present, both traditional and optimized design methods need to consider the 

above mentioned capability index. The author will put forward an optimized design 
method based on the prior knowledge. According to the method, people only need to know 
the structure of a filter and to master an intelligent optimization algorithm before finishing 
the filter's design. 

For the signal frequency of the robot joint force/ position sensors is rather low, their signal 
fits to be processed by lowpass filters. There are two kinds of filters: analog filters and 
digital filters. Here, both analog filters and digital filters are used to process the signals of 
the robot joint force/ position sensors. The configuration of the filters sees Fig. 1. 



Fig. 1. Configuration of the filters 

The output signals of the robot joint force/position sensor are analog input signals V { (V) 
of the signal processing system. After analog filtering V f (V) were converted to , and 

then were sampled and discretized into input sequence V t (ft) by A/D converter. 


2.2 Realization of the signal processing system 

(a) Realization of the analog filter 

In this research, the sensor signal is magnified by instrument magnifier AD623, and the filter 
method by double capacitors is adopted which recommended by AD623 user's manual. The 
schematic of the analog filter is shown in Fig.2. 


Vout V'out 



Cl 

+ C2 


= 0.1uF> 

— .10uF 


V V 

GND GND 


Fig. 2. Schematic of the analog filter 


(b) Realization of the digital filter 

Generally, the system function of N-order digital filter is 


Optimization Design Method of II R Digital Filters for Robot Force Position Sensors 


83 


H ( z )=— 
v ’ 1+ 


b 0 + b { z 1 +b 2 z 2 -4 f b M z M 

1 + fljZ + u 2 z + • • • + ClpjZ N 


(1) 


Translating it to difference equation 


y(n) = b 0 x(n) + b 1 x(n-l) + b 2 x(n-2)-t b a M x(n - M) 

~a 1 y(n-l)-a 2 y(n-2) a N y(n-N) 


(2) 


Then the digital filter can be realized via Eq. (2). 

3. Optimization model of the MR digital filter of robot joint force/position 
sensor 

The system of HR digital filter can be shown as Fig. 3 



HR digital filter ►Filtered signals 


Fig. 3. Schematic diagram of the HR digital filter 

Suppose that the system function of N-order HR digital filter is 



(3) 


If Eq. (3) is adopted to design HR digital filter, the number of parameter required optimize is 


M + TV + 1 , and it is difficult to choose the value range of every parameter. Generally, the 
system function of HR digital filter is expressed as 



(4) 


Both Butterworth filter and Chebyshev filter can be denoted as the cascade structural form 
with second-order unit shown as Eq. (4). When HR digital filter is denoted by this structural 
form, the sensitivity of its frequency response to coefficient change is lower. And it is 
convenient to confirm the value range of every parameter with this structure form. 

For robot force/position sensors, its measurement signal is low frequency, generally below 
10Hz. And the disturbance is white noise mostly. If power supply is mains supply, the 
disturbance of 50Hz power frequency would exist. Generally, the analog lowpass filter is 
used to deal with these kinds of signals. However, it is very difficult to filtering the 
disturbance of 50Hz power frequency and the low-frequency white noise. If the digital filter 
is adopted and its cutoff frequency is set rather low, the filter can remove the disturbance of 
50Hz power frequency and white noise mostly. From practical experience, it was known 
that the satisfying effect can be obtained when adopting a second-order lowpass. 

The system function of the second-order Butterworth filter can be simplified as 
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H(z) = A 


1 + 2 z 1 + z 2 

1 - cz” 1 + dz ~ 2 


(5) 


HR digital filter as a system, the ideal signal after filtered could repeat the input signal 
perfectly, and with definite system delay. Then the similarity of the actual and ideal output 
signal would be considered as the evaluation index of filter performance. If the maximum 
frequency of robot force/ position sensor signal is f, and the input signal is simulated with 
sine function, the function of input signal can be denoted as 

x(n) = c s •sin(2n/-n-T s ) + c dis • sin(2n- 50 -ft -T s ) + c w -randn(l) (6) 


where c s is the coefficient of sensor signal, Cdi S is the coefficient of disturbance of 50Hz power 
frequency, c w is the coefficient of white noise, T s is the sample time of system. randn( ) is the 
normal distribution random number which represents white noise. 

The ideal output after digital filter is 

y d(») = c s - sin ( 2 n /-»- r s) (7) 

The actual output of filter is 

y(n ) = A- x(n) + 2 A- x(n-l) + A- x(n-2)-c- y(n-\)- d ■ y(n -2) (8) 


Suppose that the delay time is T, the sample size is n, the mean square error (MSE) of HR 
digital filter at every sample points can be shown as 

£= ^rr2/y( T+ 0-yd(0) 2 

The mathematical model of optimization design of HR digital filter is 

fminE = /(X) 

js.f.XeS = {x|g ; .(X)<0;' = l,---,m} 
where X is optimization variable, Jf = {A,C,dj , g ; (V) is restriction function, S 
restriction region, E optimization aim function, E = f (V) . 


(9) 


(10) 


4. Optimize the parameters of the HR digital filter using the particle swarm 
optimization algorithm 

4.1 Introduction of particle swarm optimization 

Particle Swarm Optimization (PSO) algorithm is a new global optimization evolutionary 
algorithm invented by Doctor Eberhart and Doctor Kennedy. This algorithm simulates 
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preying of bird. It is outstanding to solve nonlinear optimization problem, and it is simply to 
realize this algorithm. So it has become an important optimization tool. 

In PS O algorithm, the position of particle represents the possible solution. The superiority- 
inferiority of particle is measured by particle fitness. Firstly, a flock of random particles are 
initialized. Then the optimal solution is found out via multiple iterating. During every 
iteration, particle is updated by tracking the two optimal solutions which one optimal 
solution is found by this particle, namely, individual optimal solution, the other is found by 
the whole particle swarm presently, namely, global optimal solution. 

After found two above-mentioned extremums, particle's velocity and position are updated 
according to two equations as follow. 

vr =w-V* + c 1 -r 1 \p\- x f ) + c 2 • r 2 • [p\ - xf ) (11) 

x^=x\ + V^ ( 12 ) 

where k is generation number, i is serial number of the particle, V is velocity of the particle, x 
is current position of particle, p is vector form of pBest, g is vector form of gBest, n and ri are 
random numbers from 0 to 1 , c\ and C2 are learning factors, generally, ci=C2=2, w is weighting 
factor, its value from 0.1 to 0.9. 


4.2 Improved algorithm 

(a) During the basic PSO algorithm search in solution space, the particle would oscillate 
round global optimal solution at later period. To solve this problem, the improvement can 
be done as follows: with iteration going, to let velocity update the weighting factor, the 
weighting factor decreases from maximum w m ax to minimum w m m linearly, namely 


-•G 


(13) 


where G is generation number of current iteration, G ma x is total generation number of 
iteration. 

(b) To make the PSO algorithm search in solution space and to ensure convergence rate, the 
position space and velocity space of particle need to be limited. Eq. (5) will to be 


X = X „ X , X, • , X, , Xa • , Xo 
r V 1 mm 9 1 max 9 2 mm 9 2 max 9 3 mm 9 3 max 

min 9 max 9 ^2 min 9 ^2 max 9 ^3 min 9 ^3 max 


) 

) 


(14) 

(15) 


4.3 Parameters coding and the choice of the fitness function 

To solve X using PSO algorithm, the optimization variable X should be coded, to become 
particle of PSO algorithm. According to characteristics of PSO algorithm, parameters can be 
denoted with real number. To every particle of filter shown in Eq. (5), if the current position 

of particle is denoted with X = {A,C,d} , and the velocity is denoted with 
V — ? V 2 , V 3 } , then the coding structure would be adopted as follows: 
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A, c, d 


V lr V 2/ V 3 


Position of the particles Velocity of the particles 
PSO algorithm confirms the superiority-inferiority of particle's current position via fitness. 
So the fitness function must be chosen according to the practical demand. Here Eq. (9) is 
chosen as the fitness function of HR digital filter design. Evidently, the less the value of E is, 
the less the mean square error of filter parameter X corresponding this particle is. Then this 
particle is corresponding to the better coefficient of filter. When the algorithm finished, the 
particle with the minimum fitness during the whole running period is the optimal solution 
obtained by this algorithm, namely, filter parameter. 


4.4 Optimization steps 

(a) Set parameters of PSO algorithm, including population size S pop , dimension Sdim, 
weighting factor w, position space x r , velocity space v r . 

(b) Set parameters of HR digital filter, such as c S/ Cdi S , c w , T S/ etc. 

(c) Initialize the particle swarm, to randomly initialize every particle's position and velocity 
in parameter space. 

(d) Solve the system delay time, and to calculate the particle's fitness according to Eq. (9). 

(e) Initialize the current particle's position as individual extremum /?Best, and the position of 
particle with minimum fitness among all individual extremum as gBest- 

(f) Update the particle's position and velocity according to Eq. (14) and Eq. (15). 

(g) Solve the system delay time, and to calculate the particle's fitness again. 

(h) Judge whether to update the particle's individual extremum psest and the global 
extremum gBest of particle swarm. 

(i) Repeat step (f) to (h), till meeting precision demand or reaching iteration times pre-set. 

(j) Output gBest, and to obtain the coefficient of HR digital filter. 


5. Results and analysis 

To prove the validity of optimization design method of HR digital filter for robot 
force/ position sensor presented in this paper, an optimization program was completed in 
the Matlab 6.5 circumstance and several simulation experiments were made. In these 
simulation experiments. Parameters of the PSO algorithm were set as: swarm size S pop = 50, 
parameter dimension Sdim = 3, minimum weighting factor zc min = 0.1, maximum weighting 
factor w max =0.9. According to former design experience, the value range of every filter 
coefficient are as: x r = [0.00001, 0.01, 1, 2, 0.00001, 1], maximum velocity v T = [-0.005, 0.005, - 
0.5, 0.5, -0.5, 0.5], sample time T s = 0.02 s, maximum frequency of sensor signal /= 0.05 Hz, 
signal amplitude c s = 1.5, disturbance amplitude of 50 Hz power frequency Cdi S = 0.1, 
amplitude coefficient of white noise c w = 0.15, simulation time t = 5 0 s. 

Fig. 4 is the variety course of every generation's fitness of PSO algorithm. 
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Fig. 4. Curves of unfiltered and filtered signals of the robot force/ position sensors 


Fig. 5 a) and b) are the signal curve of robot force/ position sensor before and after filtering 
respectively. The last parameters of HR digital filter are set as A = 0.0064, c = 1.5186, d = 
0.5439, delay time T = 0.36 s, MSE E = 0.0091. From Fig. 4 b), we can know that the filter 
effect is perfect for sensor signal of /= 0.05 Hz. 



(a) (b) 

Fig. 5. Curves of unfiltered and filtered signals of the robot force/ position sensors, (a) is 
curve of unfiltered signals, (b) is curve of filtered signals. 


The frequency response optimized is 


fi(e T >') = 0.0064 x 


1 + 2e~ 0<a ’°’ + (e~ 0 ' 02 '"") 2 

1 - 1.5186e~°' 02;< " + 0.5439(e _0 ' 02;< ") 2 


(1 

6) 


From Eq. (16), the amplitude-frequency characteristics curve of filter designed in this 
research can be obtained, as shown in Fig. 6. 



°0 OS I 1$ 3 25 3 35 4 4$ 

/(Hz) 

Fig. 6. Magnitude-frequency characteristics curve of the filter 
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From Fig. 6, it can be seen that signal is attenuated under 0.1, when its frequency decreased 
to 5 Hz. So the filter can remove the disturbance of 50 Hz power frequency and white noise 
mostly. 


6. Conclusion 

Aimed at the design of HR digital filters of robot force/ position sensors, a design method is 
put forward. Its optimization principle is the minimum MSE between ideal and actual 
output signal at time-domain. And the mathematics model aiming at second-order 
Butterworth lowpass filter was set up. This method needn't understand the complicated 
design theory and method for digital filter and the characteristic of filter, such as passband 
frequency, cutoff frequency, passband attenuation, ripple, etc. This method only requires the 
understanding of the structure characteristic of filter and the maximum frequency of sensor 
signal. Thus the parameters of the filter can be optimized in a suitable intelligent 
optimization method. An optimization program of the PSO algorithm was developed in the 
Matlab circumstance. The result of simulation experiment proves the validity of this 
method, and to be strongly practicable. 
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Visual Analysis of Robot and Animal Colonies 
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1. Introduction 

Micro-robotics calls for the development of tracking systems in order to study the 
movement of each micro-robot in a colony to answer questions about what they are doing 
and where and when they act (see Fig. 1). Moreover, micro-robots can be used to emulate 
social insect behaviour (Camazine et al., 2001) and study them through tracking experiments 
involving several miniature robots on a desktop table. Thus, principles of self-organization 
in these colonies, which were studied so far by analysis of a tremendous amount of insect 
trajectories and manual event counting, are now better understood by biologists thanks to 
robotics research. 



Fig. 1. Analogy between social insects (left) and a micro-robot colony (right) 

Although this approach is only recently increasing its popularity, computer vision systems 
for tracking moving targets are widely used in many applications such as smart 
surveillance, virtual reality, advanced user interfaces, motion analysis and model-based 
image coding (Gavrila, 1999). Surveillance systems seek to automatically identify people, 
objects or events of interest in different kind of environments (Russ, 1998; Toyama et al. 
1999; Haritaoglu et al., 2000; Radke et al., 2005). However, this problem is not easy to solve. 
First of all, it is not viable to tag each colony member under study. On the one hand, the tag 
selection process can be difficult since tags must be very small in some cases and, therefore, 
it might not be possible to detect them in an image. Furthermore, tags can become 
ambiguous when a swarm is composed of many individuals. On the other hand, tagging can 
alter individual behaviour. So, an application for tracking unmarked object has been 
developed. A new problem arises: how to identify the same object in two consecutive frames. 
SwisTrack (Correll et al., 2006) is a previous work following this approach which we try to 
improve. It is a platform-independent, easy-to-use and robust tracking software developed 
to study robot swarms and behavioural biology. It is part of the European project LEURRE 
(http:/ /leurre.ulb.ac.be.2006) focused on building and controlling mixed societies composed 
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of animals and artificial embedded agents. In preliminary case studies towards this aim 
(Caprari et al., 2005), an Insbot team has been introduced into a swarm of cockroaches and 
allowed for modification of the natural behaviour of the swarm. We will show that our 
application achieves robust performance in object identification and tracking without the 
need of a strong intervention by the user. 

As the structure of the designed method (see Fig. 2), this paper is organized as follows: the 
visual segmentation and detection of objects are described in Sec. 2 and 3. In Sec. 4 we 
outline the tracking problem and its solution. Experimental results are given in Sec. 5 and 
discussed in Sec. 6. 



Fig. 2. Flowchart of the whole designed method 


2. Object Segmentation 

A common element of any surveillance systems is a module that is capable of identifying the 
set of pixels which represents all the individuals under study in each captured image. There 
are several techniques to carry out this task. For example, the background modeling 
approach (Toyama et al., 1999; Haritaoglu et al., 2000) allows to model dynamic factors such 
as blinking of screens, shadows, mirror images on the glass windows or small variations in 
illumination due to flickering of light sources. Pixels are classified as background or 
foreground depending on the fitting of their values with the built model. As a drawback, 
this method is not capable of adapting to sudden illumination changes, and we ruled it out 
with the aim of developing a more robust surveillance application in the presence of 
variation of lighting conditions. On the other hand, most of the alternative techniques are 
developed by gray-level image processing so that if available images are in color, it is 
necessary to convert them to gray-level. As our system has color images as input, the binary 
image resulting from the segmentation process can be obtained from a combination of three 
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binary images (each color channel generates a gray-level image which is segmented by the 
selected method obtaining a different binary image), or from a gray-level image obtained 
directly from the color one. 

Thus, the first step is to convert the captured color image to a gray-level one. For such 
preprocessing, the Hue-Saturation-Intensity (HST) system is used since it encodes color 
information by separating an overall intensity value I from two values encoding chromaticity 
- hue H and saturation S. HSI might also provide better support for computer vision 
algorithms because it is amenable to normalization for lighting and focus on the two 
chromaticy parameters that are more associated with the intrinsic character of a surface 
rather than the lightning source. Thus, the resulting gray-level image is only built from the 
intensity value. Derivation of HSI coordinates from RGB coordinates is a common process in 
computer vision (Shapiro and Stockman, 2001). 

Once the gray-scale image is available, a segmentation process has to be applied on it. 
Although frame difference is the easiest and fastest method to detect moving objects in an 
image, it fails when the objects are steady. This problem could be solved by taking a 
reference image with no objects and subtracting it from the new ones. Nevertheless, a little 
illumination change might make the whole process fail, thus, an alternative algorithm 
should be chosen. In our case, the corresponding binary image is obtained from an input 
gray-scale by thresholding operations as in Swistrack. This technique defines a range of 
brightness values in the original image: pixel value greater than a threshold (or lower, 
depending on its definition) belongs to the foreground and the rest of pixels are classified as 
background. The drawback of this method is the correct determination of the threshold. In 
Swistrack two different kinds of reference images are used, depending on the mode in which 
the system is operating: 

• Static background: the background image is captured at the beginning of the 
experiment and is not updated at any time; 

• Running average: the reference image is built as the running average of all video frames 
processed until that moment 

The threshold is fixed in both cases and represents the minimum difference required to 
classify a pixel as foreground. It is important to note that the fixed threshold is sensitive to 
changes in lighting conditions, especially in the first operation mode in which the reference 
image is not updated during all the experiment. Although the running average is more 
capable of dealing with illumination changes, it might consider objects that stop moving for 
a long period of time as part of the background, without detecting their presence in the 
scene. We have implemented a method for automatically calculating the threshold based on 
histogram properties by updating its value in each frame, in order to adapt it to variations of 
the lighting conditions. This provides an advantage over the Swistrack method, as significant 
intensity differences between the objects to be tracked and the background are not 
necessary. 

After the threshold setting, two consecutive morphological operations are applied on the 
binary image resulting from the segmentation process. These steps are required to erase 
isolated points or lines caused by different dynamic factors such as, for example, changes 
induced by camera motion, sensor noise, non-uniform attenuation, blinking of lights or 
atmospheric absorption. A 3x3 erode filter is used to delete these artefacts and then a 3x3 
expand filter is applied to recover the foreground region. A result of the whole process can 
be observed in Figure 3, where a group of micro-robots is seen from above. As it can be seen. 
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although the lighting conditions are not good in some places, the designed application is 
capable of detecting all visible micro-robots in the image. However, the pixels due to light 
reflexions on the arena are not removed from the image. This issue will be solved by means 
of the implemented labeling method described in next section. 



Fig. 3. An image captured by the system camera (left), the binary image obtained by the 
segmentation process (center) and the resulting binary image after applying two 
morphological operations (right) 


3. Object Identification 

The aim of this stage is to obtain a labeled image in which each label identifies one colony 
member. This can be a difficult task when objects are touching one another, because an 
identified blob contains all objects with, at least, one point in common. The method 
implemented to achieve the above goal can be divided into three steps: 

1. Labeling of the identified blobs in the input binary image. A row-by-row labeling 
method (Shapiro and Stockman, 2001) is used to reduce the computational cost of the 
whole process. The binary image is scanned twice: the first time to tag each 
foreground pixel based on the labels of its neighbors and to establish equivalences 
between different labels; the second, crossed scan, will unify tags which belong to the 
same blob 

2. Classification of the labeled blobs as targets to be tracked or as bad-segmented pixels, 
and detection of collisions inside a blob identified as a tracking target. All targets are 
assumed as not touching in the first captured image. The number of targets to be 
tracked is specified by the user, and the application calculates the minimal and 
maximal size allowed from the first captured image. This knowledge, together with 
the number of the detected blobs in each frame, allows to define a series of criteria to 
determine when a blob represents more than one object, and to reject all blobs that do 
not identify target objects but are instead the result of a bad segmentation, as it occurs 
with the set of noisy pixels in the example of Fig. 3 

3. Segmentation of blobs that represent groups of more than one object to be tracked. This 
step is explained in more detail in the next section 

As seen above, it is possible that the same blob identifies more than one object to be tracked. 
Thus, it is important to detect these situations and split up the blob in the corresponding 
objects. There are two different tasks to be performed: determining the number of touching 
objects inside the same blob and splitting them up. 

The first goal is achieved through several criteria which are relationships between blobs and 
object sizes. They can be easily set up by the application, assuming that no target objects are 
touching in the first captured image. Thus, our application calculates the maximum and 
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minimum dimensions of the visual objects from the first captured scene and the criteria 
remain set. This step is important because perceived object size can vary with the distance 
between the arena and the camera and if an object size in pixels is directly related to its real 
size, as in Swistrack, a calibration error can be introduced in all calculations. Due to this 
error, the application might fail the tracking process. So, the only parameter our application 
needs, which is requested to the user, is the number of objects to be tracked. 

The next step is to split up the different objects which compose one blob. As it is difficult to 
identify several objects at the same time, we have studied three different, possible situations 
assuming that only two objects are touching. Thus, our application split up any complex blob 
in two different parts: one target object and another blob which can be again composed of 
more than one micro-robot. If the new blob represents several micro-robots, the split-up 
process is recursively applied until the obtained blob is composed of only one target object. 
The three different cases that have been studied are the following: 

1. two objects touching only in one point; 

2. two objects sharing one side, that is, they are horizontally or vertically aligned; 

3. two objects touching in several points which do not correspond to their sides 

In the first case, a contour method is used. A chain code is calculated by considering that the 
contact pixel will be visited twice. The method takes into account the contour irregularities 
due to the segmentation process and it applies different criteria to determine the correct 
contact point as shown in Fig. 4. 



The second case, when two objects share one side, has been considered because of the micro- 
robots used in our experiments, the Alice2002 (Caprari and Siegwart, 2005) (see Fig. 5). These 
micro-robots can be seen from above as boxes and two Alice2002 can share one side in any 
moment. It is important to note that this case does not apply to micro-robots that do not 
have a shape in which a side can be shared. A method based on dimension criteria is used to 
determine the common side, and it estimates their splitting line as shown in Fig. 4. 



back 


(a) A micro-robot Alice 2002 (b) Infrared proximity sensors 

Fig. 5. Micro-robot Alice 2002 

Finally, the most general and complex case is when an object and another blob compose a 
bigger blob, and the contact between them is through several pixels which do not correspond 
to an object side. Therefore, the designed method is based on holes inside blobs. Again, this 
might be the result of a bad segmentation. For this reason, a set of criteria were defined for 
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detecting when a hole is due to bad segmentation and when it is a hole between two 
different objects. As it can be seen in Fig. 4, the designed method provides successful results. 


4. Tracking Micro-robots 

Once each micro-robot is identified as an object, the next step is to match each object with 
one of those detected in the previous frame, in order to obtain its trajectory. 

Object position can be calculated as the geometrical center of gravity of object contour 
(Correll et al., 2006), but we decided to calculate it as the geometrical center of gravity of a 
blob corresponding to an object. Since the correspondence between an object and a blob has 
been obtained in the previous step, this procedure is easier and faster. The issue now is how 
to associate each center of gravity with its corresponding object between the new ones. The 
chosen solution, the nearest neighbor method, is an easy one, but its implementation raises 
several issues (Correll et al., 2006). The different way of dealing with them will determine 
the successfulness of the application. 

The first challenge to solve is the case of an object that is closer to the previous position than 
the real one (situation already described in a previous work (Correll et al., 2006)). A 
quadratic assignment problem for minimizing the sum over the distance of all assignments 
is used in Swistrack, but this does not constitute an optimum solution since it fails in some 
cases. On the contrary, a solution focused on the previous movements of objects is presented 
here, that is, the matching algorithm associates the information on the nearest neighbor with 
that regarding the previous movement direction. 

It might also be that an object disappears from the scene (Correll et al., 2006), but if the 
application does not detect all elements specified by the user, then it will fail when this 
situation occurs, because objects enter or leave the arena, and the system will repeatedly 
capture another image until it finds all the objects. Another situation (again presented in 
(Correll et al., 2006)) is when two shared trajectories are divided at the wrong time, but this 
is not possible in our application because blobs are divided into their corresponding objects 
even though they are touching, as previously explained. Finally, an additional skill of our 
application is the ability of avoiding un-associated contours. 

Overall, our application is capable of solving the four different situations expounded in 
(Correll et al., 2006) and does not need the help of the user for correcting wrong trajectories. 
It is also important to note that the modified nearest neighbor technique is only applied on a 
small region of the image, not on the whole image, in order to make the application faster. 
Again, it is not possible to fix the search area size because a delay can be produced during 
capture process. Thus, our application calculates the search area dimension based on 
capture time between frames and the maximum velocity of tracking objects, which will be 
requested to the user. This reduces the computational time and guarantees the success of the 
matching process. 


5. Experiments 

We first provide a brief overview of the robotic platform used, followed by the experimental 
results. 

The experimental setup is depicted in Figure 6. A camera is pointing downwards to the 
desktop where the micro-robots are working. Our application operates with monocular 
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color video images (see Fig. 6). The distance between the camera and the arena can vary 
from one setup to another and the system calculates the parameter values needed for 
obtaining the different criteria. It is important to note that the background of the micro-robot 
workspace is black in all our experiments, while the most common solution is to use a white 
background. 



Fig. 6. Experimental setup (left) and a 320x240 captured image by the camera (right) 

The objects to be tracked are Alice 2002 robots, as mentioned above. They are extremely 
sensitive to external forces and can be very easily damaged if not handled with care. Among 
their features, we can mention: 

• tiny dimensions (22 mm (width) x 21 mm (length) x 20 mm (height)) 

• small weight (aproximately 11 g) 

• infrared proximity sensors (front, right and left) to avoid obstacles 

• low consumption (12 - 17 mW) 

• autonomy (up to 10 hours thanks to its Ni-MH rechargeable battery (Varta 3/V40H)) 

• velocity of 40 mm/ s 

Finally, a Graphical User Interface (GUI) has been developed to check the performance of 
our application. It is composed of two different windows (see Figure 7(a)): on the left, the 
user can observe the images taken in real-time and, on the right, a graphical window is 
showing the different positions of the objects. Each obtained trajectory is drawn in a 
different color to help the user identify each target. As the duration of the experiments is 
unknown and the amount of points can be considerable, the application only shows the last 
seventy object positions to ease tracking of each described trajectory to the user. 

Two different experiments have been carried out. The first one is the tracking of three 
unmarked AlicelOOl (see Figure 7). Six untagged Alice2002 are studied in the second 
experiment (Fig. 8). Both experiments have been carried out at different times of the day and 
in different days to test the robustness of our application to different lighting conditions. 

As it can be seen in Fig. 7a, all micro-robots to be tracked are not touching in the first stage. 
Thus, our application can obtain the information it needs: the objects position, i.e., their 
geometrical center of gravity, and the maximum and minimum size allowed for any object. 
For clarity of representation, objects are highlighted by inscribing them in circles. 

Although there are relevant delays between the first and the second frames and between the 
second and the third ones, our application is capable of correctly tracking the objects as 
shown in Fig. 7b. 

To check the system in different situations, we have changed the moving pattern of the 
objects during the experiment. For the first 25 frames all three objects are describing a 
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clockwise circular trajectory. A different circular trajectory with smaller radius is described 
during the next 25 frames, as shown in Fig. 7d. Finally, the micro-robots are set in wall 
following mode, so they describe a straight-line trajectory until they find a wall to follow, as 
can be observed in Fig. 7e. 

The second experiment was similar. In this case, our application had to track six unmarked 
Alice2002. Again, the first captured frame (see Fig. 8a) reveals that the experiment starts 
without collisions between objects. The trajectories described in this experiment are, first, 
the smallest-radius circular trajectory; then, four of the micro-robots change their trajectory 
describing a circular trajectory with a larger radius and the two remaining ones go straight 
ahead looking for a wall to follow. Now, two more micro-robots change to the wall 
following mode. It is important to note that our application is capable of detecting objects 
only partially visible as shown in Fig. 8d. 



. \ T ... i c , . _ , (c) Circular (e) Following wall 

(a) Initial State (b) After 5 frames . trajectory with a v ’ , ° 

v ' v ' trajectory ' ' mode 

' J different radius 

Fig. 7. Three Alice 2002 tracking experiment: captured image (up) and trajectory image 
(down) 
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Fig. 8. Six Alice 2002 tracking experiment 


In addition, examples of the use of our split-up method are finally shown in Fig. 9. A one- 
shot video segment of 10 minute duration is available at 
http://www.robot.uji.es/lab/ plone/ Members/ emartine 
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Fig. 9. Collision detection 


To conclude this section, a graph is presented (see Fig. 10) which compares the trajectory 
followed by a member of the studied colony in a multiple-target experiment versus the 
trajectory obtained by the developed software. As it can be observed, there is no mismatch 
in data-association thanks to the implemented method to split up blobs when several 
members are in touch. 


Trajectory followed by an alice vs the one obtained by the developed application 

I Trajectoiy followed I 

• i | T rajectory calculated | 



Trajectory Error 



Fig. 10. Trajectory followed by a single target vs individual trajectory obtained by the 
implemented software (above) and error (below) 
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6. Conclusions 

We have presented a tracking application to study micro-robots or social insect cooperative 
behavior without the risk of conditioning the results by tagging them. Our system has been 
compared with previous ones, and namely with Swis track, an application intended to control 
mixed societies. Although this previous study had the same goal, the authors deal with the 
tracking problem in a different way. The given results have shown the robustness of our 
application with regard to lighting conditions. Also, no special illumination is required and 
performances do not depend on the surrounding objects, as for example it occurs in 
Swistrack. 

Our designed method also solves situations in which there are several objects touching one 
another and it can match an object position in one frame with its position in the next frame. 
It is also capable of detecting objects even though their velocity is very slow or if they do not 
move, a case typically difficult for similar methods. 

As a further achievement, our application only requires two parameters from the user: the 
number of target objects and their maximum speed. No thresholds need to be set manually. 
Overall, we have designed an application transparent to the user who does not need to 
know the implementation details to work with it. 

So far, our application has only been tested with homogeneous robotic societies. As further 
research, we plan to test it with mixed societies. 
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